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■* 

O&JECTI VES  OF  THE  RESEARCH  PROGRAM 

The  main  emphasis  of  the  research  program  is  placed  on  the  dynamical 
performance  of  truss-type  space  structures,  with  the  following  broad  objectives: 

(1)  Develop  a  basic  theory  capable  of  perdicting  accurately  the  response 
of  a  truss-type  space  structure  under  arbitrary  loading,  without  making  the 
usual  sweeping  approximation  in  which  actual  truss  networks  are  replaced  by 
Timoshenko  beam  elements  or  other  continuum  elements  in  the  analysis. 

(2)  Evaluate  the  effects  of  random  variations  in  material  and  geometrical 
properties  in  the  truss-network  construction. 

(3)  Identify  nonlinearities  in  aerospace  strucutural  systems  and  develop 
new  exact  and  approximate  solution  techniques  for  such  problems. 

DETERMINISTIC  ANALYSIS  OF  TRUSS-TYPE  STRUCTHRA1  NEIMQRJl 

Truss-type  networks  are  often  used  in  space-structure  designs.  Their 
configurations  are  generally  complex,  and  a  large  number  of  degrees  of  freedom 
is  required  to  adequately  describe  their  responses  to  dynamic  excitations.  The 
application  of  an  all-encompassing  finite  element  analysis  can  be  expensive  in 
computer  time.  In  many  cases,  it  may  also  lead  to  inaccurate  results  since  some 
of  the  natural  frequencies  may  be  c.v  close  to  each  other.  One  way  to 
alleviate  this  predicament  is  to  subst,.jte  the  original  truss-type  network  by 
an  equivalent  network  of  Timoshenko-beam  elements  (e.g.  references  [1,2]). 
However,  much  detailed  information  is  lost  in  the  process  which  may  be  important 
for  certain  missions  of  the  structure  or  important  in  the  design  of  control 
systems  for  the  structure. 

Recognizing  the  fact  that  a  typical  truss  network  of  a  space  structure 
consists  only  of  a  few  sizes  of  truss  units  which  are  connected  end-to-end  to 
form  piece-wise  periodic  arrays,  recent  advances  in  the  field  of  periodic 
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structures  (e.g.  references  [4-12])  can  be  utilized  to  advantage.  The  present 
analysis  is  bas^d  on  the  concept  of  wave  propagation  in  periodic  structures  [4], 
modified  to  account  for  wave  scattering  when  a  periodic  pattern  is  distrupted. 
The  finite  element  or  some  other  conventional  formulations  may  be  used  to  model 
a  single  truss  unit,  transfer  matrices  are  used  to  relate  the  state  vectors  of 
structural  response  from  unit  to  unit,  and  state  vectors  are  converted  to  wave 
vectors  in  order  that  computation  can  follow  the  direction  of  wave  propagation 
to  avoid  numerical  instability  in  the  computation.  The  general  approach  is 
similar  to  that  adopted  by  von  Flotow  [13]  who  substituted  the  truss  network  by 
Timoshenko-beam  network,  and  by  Signorell  and  von  Flotow  [14]  who  were  concerned 
with  wave-passage  frequency  bands,  wave  modes,  and  power  flow  in  a  purely 
periodic  truss  beam.  A  combined  finite  element  and  transfer  matrix  approach  was 
first  proposed  by  McDaniel  and  Eversole  [15].  The  retention  of  truss  unit  and 
the  inclusion  of  piece-wise  periodic  patterns  and  structural  damping  in  the 
analysis  do  add  considerable  complications;  however,  the  rewards  are  more 
accurate  results,  and  broader  applicability. 

1.  MULTI-CHANNEL  WAVE  GUIDES 

Any  identifiable  structural  element  may  be  viewed  as  a  waveguide,  in  the 
sense  that  disturbances  at  one  location  can  be  transmitted  to  another.  Consider 
a  typical  truss  unit  in  a  space  structure,  which  has  an  exterior  shape  of  a 
right  parallelepiped,  with  eight  exterior  corners.  If  the  generalized 
displacements  and  the  associated  generalized  forces  at  the  four  corners  on  one 
end  of  the  unit  are  specified,  then  the  generalized  displacements  and  the 
associated  generalized  forces  at  the  four  corners  on  the  other  end  can  be 
determined.  Since  each  corner  is  capable  of  six  degrees  of  freedom,  three  in 
translation  and  three  in  rotation,  such  a  unit  is  capable  of  transmitting  24 
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types  of  wave  motions  if  the  structural  members  in  the  unit  are  rigidly 
connected  at  the  corners.  In  the  terminology  of  wave  mechanics,  it  is  a 
24-channel  waveguide.  However,  if  the  structural  members  are  pinned  at  the 
corners,  then  the  rotational  degrees  of  freedom  cannot  be  transmitted,  and  the 
unit  reduces  to  a  12-channel  waveguide. 

To  be  more  specific,  write  an  equation  in  the  frequency  domain  relating  the 
state  vectors  on  the  two  ends  of  a  truss  unit,  say  unit  n  (or  cell  n),  as 
fol lows 


rwn+l.  rwn, 

{  }  =  [T (n)]  {  } 

fn+l  fn 


(1) 


where  each  w  is  a  p-dimensional  vector  of  the  complex  amplitudes  of  generalized 
displacements,  each  f  is  a  p-dimensional  vector  of  the  complex  amplitudes  of  the 
generalized  forces,  and  T  is  a  2px2p  transfer  matrix  associated  with  the  nth 
cell.  For  a  typical  truss  cell,  p=24  or  p=12  depending  on  whether  the  truss 
members  are  clamped  or  pinned  at  the  corners.  Using  the  transformation 


wn  un(n.n+l) 

{  }  «  [D(n)]  {  } 

fn  pn(n+l,n) 


(2a) 


rwn+l.  Mn+l(n,n+l) 

{  }  «  [D(n)]  {  } 

fn+1  Mn+l(n+l»n) 

where  columns  in  [D(n)]  are  the  eigenvectors  of  [T(n)]  arranged 
ascending  order  of  the  magnitudes  of  the  eigenvalues,  equation 


(2b) 

according  to  an 
(1)  is  replaced 


by  [3] 


^n+l(n»n+l)  ^ 

pn+l(n+l,n) 


A  0  '  ^n(n,n+l) 
0  A-1]  yn(n+l,n) 


(3) 


in  which  A  is  a  diagonal  matrix  of  the  first  p  eigenvalues.  When  deriving 
equation  (3)  use  is  made  of  the  fact  that  the  eigenvalues  of  a  transfer  matrix 
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are  reciprocal  pairs.  Since  damping  normally  exists  in  a  physical  system,  the 
magnitude  of  each  eigenvalue  in  A  is  smaller  than  one,  and  each  in  A"1  is 
greater  than  one.  This  implies  that  yn(n,n+l)  in  the  first  row  of  equation  (3) 
represents  a  wave  vector  transmitted  through  cell  n,  and  emerges  as  yn+l(n,n+l). 
The  eigenvalues  in  A  account  for  both  reduction  of  magnitudes  and  change  of 
phases  in  the  transmission  process.  In  the  second  row  of  equation  (3),  it  is 
the  wave  vector  yn+i(n+l,n)  which  is  transmitted  through  cell  n  and  emerges  to 
become  yn(n+l,n).  The  choice  of  the  subscript  and  the  order  of  the  arguments  in 
denoting  a  wave  vector  y  can  now  be  explained.  For  example,  yn(n,n+l) 
represents  a  wave  vector  at  location  n  which  is  propagated  from  n  to  n+1,  etc. 
The  cell  plays  the  role  of  a  p-channel  waveguide. 

Obviously,  two  p-channel  waveguides  connected  end-to-end  is  also  a  p- 
channel  waveguide.  Consider,  for  example,  two  connected  cells,  and  write  the 
transfer  relationship  for  the  state  vectors  as  follows: 


_wn+2,  wn 

{  )  ■  [T(n+l)][T(n)]  {  J 

fn+2  fn 

It  follows  upon  transformation  into  wave  vectors 

yn+2(n»n+2)1  yn(n,n+2) 

{  I  =  [Q(n+2,n)]{  } 

yn+2(n+2*n)  yn(n+2,n) 


(4) 


(5) 


where 

[Q(n+2,n)]  =  [D(n+l)]-l[T(n+l)][T(n)][D(n)]  (6) 

Now  [0(n+2,n)j  is  generally  not  a  diagonal  matrix,  unless  [T(n+l)]=[T(n)j, 
namely,  unless  cells  n  and  n+1  are  identical.  It  implies  that  yn+2(n,n+2)  does 
not  result  completely  from  the  transmitted  portion  of  yn(n,n+2). 


2.  WAVE  SCATTERING  MATRICES 

For  an  arbitrary  array  of  connected  cells,  equation  (5)  is  generalized  to 
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Vm(m.n) 

V*m(n,m) 


Un(m.n) 

}  -  CO(n,m)]  {  } 

yn(n,m) 


(7) 


These  wave  vectors  are  shown  schematically  in  figure  1.  With  respect  to  this 
waveguide,  ^(m.n)  and  yn(n,m)  may  be  treated  as  incoming  wave  vectors  whereas 
ym(n,m)  and  yn(m,n)  treated  as  outgoing  wave  vectors.  These  two  types  of  wave 
vectors  can  easily  be  recognized  from  their  notations.  A  wave  vector  is 
incoming  if  the  subscript  and  the  first  argument  are  the  same;  it  is  outgoing  if 
the  subscript  and  the  second  argument  are  the  same. 

Recast  equation  (7)  in  a  form  indicative  of  an  input-output  relationship  as 
f o 1 1 ows : 


Vm(n.m) 

yn(m,n) 


,VmOn,n) 

=  [S(n,m)]  {  } 

yn(n,m) 


where 


[S(n,m)]  = 


R(m,n) 
T (m,n) 
It  can  be  shown  that 


T(n,m) 

R(n,m) 


(9) 


R(m,n)  =  -  [Qjjtn.m)]"1  Q21(n,m)  (10a) 

T(m,n)  =  Qn(n,m)  -  Ol2(n,m)[Q22(n,m)]-1Q21(n,m)  (10b) 

R(n,m)  ■  Ql2(n,m)  [Q22(n,m)]~l  (10c) 

T(n,m)  =  [022 (n,m)]_1  (lOd) 

where  Qij(n,m)  are  sub-matrices  of  Q(n,m)  in  equation  (7).  The  doubly 


argumented  T  and  R  matrices  in  equation  (9)  are  called  transmission  and 
reflection  matrices,  respectively.  The  first  argument  in  the  parentheses 
denoting  these  matrices  indicates  the  location  of  an  incoming  wave  vector  to  be 
transmitted  or  reflected.  For  example,  R(m,n)  indicates  a  reflection  response 
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to  an  incoming  wave  vector  at  m,  and  T(n,m)  indicates  a  transmission  response  to 
an  incoming  wave  vector  at  n.  Moreover,  the  two  arguments  for  example,  (m,n) 
refer  to  a  waveguide  extending  between  m  and  n,  and  their  order  indicates  that 
the  incoming  wave  is  in  the  direction  of  m  to  n.  Matrix  S  in  equation  (8), 
composed  of  these  transmission  and  reflection  sub-matrices,  is  called  the  wave 
scattering  matrix. 

It  is  of  interest  to  point  out  that  an  array  of  cells  may  contain  several 
sections,  each  of  which  is  composed  of  one  type  of  periodic  units.  Thus,  the 
entire  array  is  a  piece-wise  periodic  structure.  The  scattering  matrix  for  each 
section  is  easy  to  construct  since  their  transmission  matrices  are  of  the  form 
of  AJ  where  j  is  the  number  of  the  periodic  units  in  the  section,  and  their 
reflection  matrices  are  null.  The  scattering  matrix  for  the  interface  between 
two  different  periodic  sections  can  be  determined  from  a  wave-vector  transfer 
matrix  0  of  the  form 

0  =  CD(j+l)1-l[D(J)]  (11) 

according  to  the  rule  given  by  equations  (lOa-d). 

If  a  waveguide  is  terminated  at  a  boundary  where  either  the  force 
components  are  zero  (a  free  boundary)  or  the  displacement  components  are  zero  (a 
clamped  boundary),  then  such  a  boundary  can  also  be  treated  as  an  interface 
between  the  end  of  the  waveguide  and  a  fictitious  station  beyond  the  waveguide. 
For  instance,  let  n  be  a  free  boundary  such  that  Fn  =  0.  We  may  write 

fwn  r  /  xi  f^n(n»  Fn)  Mn(n,Fn)  Mn(Fn,n)l  Mn(n»Fn), 

{  }  -  CD(n)]  {  }  =  {  }  (12) 

0  Pn(Fn»n)  Nn (n ,Fn)  Nn(Fn,n)J  yn(Fn,n) 

where  Fn  indicates  a  fictitious  station  next  to  n  but  outside  the  waveguide,  and 

where  the  Mn  and  Nn  are  sub-matrices  of  D(n),  which  play  the  roles  of  converting 

the  outgoing  waves  pn(n,Fn)  and  incoming  waves  yn(Fn,n)  into  displacements  and 

forces,  respectively  [8].  From  the  second  row  of  equation  (12), 


6 


Mn(Fn»n)  =  '  CNn(Fn»n)]-1  Nn(n-Fn)Mn(n*Fn)  (13) 

Now,  the  interface  between  n  and  Fn  may  be  characterized  by  a  scattering  matrix, 

namely, 


{ 


Mn(Fn  *n) 

Mr  (n »Fn) 
h  n 


}  =  [S (Fn ,n) ]  { 


Mn(n  »Fn) 

Me  ( Fn » ^ ) 
r  n 


R(n,Fn)  T(Fn,n) 
T(n,Fn)  R(Fn,n) 


^Mn(n,Fn) 
Mr  (Fn*n) 


(14) 


However,  neither  incoming  nor  outgoing  waves  exist  beyond  Fn,  we  also  have 
MFn(Fn»n)  ■  MFn(n,Fn)  =  0.  It  follows  that  T(n,Fn)  =  0,  and 


R(n ,Fn)  =  -  [Nn(Fn,n)]-l  Nn(n,Fn)  (15) 

Matrices  T(Fn,n)  and  R(Fn,n)  can  be  arbitrary,  but  they  are  not  needed  in  our 
cal culation. 

A  clamped  boundary  where  the  displacement  components  must  vanish  can  be 
treated  similarly.  Let  Cn  be  a  fictitious  station  next  to  a  clamped  boundary  n. 
It  can  be  shown  that  at  the  boundary,  the  transmission  matrix  T(n,Cn)  *  0  and 
the  reflection  matrix 


R(n,Cn)  =  -  [Mn(Cn,n)]-l  Mn(n,Cn)  (16) 

the  concept  can  further  be  generalized  to  the  case  in  which  any  p  components  of 
a  2p  state  vector  must  vanish  at  the  boundary.  If  the  vanishing  components  are 
of  a  mixed  type,  including  some  force  components  and  some  displacement 
components,  it  will  be  necessary  to  re-arrange  the  state  vector  such  that  the 
vanishing  components  are  grouped  in  either  the  lower  half  or  the  upper  half. 
The  rows  in  the  M  and  N  matrices  must  also  be  rearranged  accordingly.  The 
reflection  matrix  then  can  be  obtained  in  a  similar  manner.  In  what  follows, 
such  boundaries  are  referred  to  as  bounded  ends. 


7 


Finally,  the  scattering  matrix  for  a  chain  of  waveguides  and  interfaces  can 
be  constructed  equentially  using  the  composition  rule  [3] 


R(i,k) 

=  R(i  ,j)  +  T ( j , i ) 

R(j,k)  [I-R(j.i)  R ( j ,k)]"l 

(17a) 

T(i,k) 

■  T(j,k)  [ I -R ( j , i 

)  R( j ,k)]-l  T ( i , j ) 

(17b) 

where  station  j  is  between  stations  i  and  k.  It  is  also  clear  from  the  above 
discussion  that  wave  reflections  in  a  piece-wise  periodic  waveguide  occur  only 
at  the  interfaces  and  boundaries  where  changes  of  periodic  patterns  take  place. 

Computationally,  the  use  of  scattering  matrices  in  the  formulation  has  a 
distinct  advantage;  it  directs  the  calculation  in  the  directions  of  wave 
propagation  and  ensures  numerical  stability.  Since  damping  exists  in  a 
structure,  the  amplitude  of  a  wave  motion  decreases  as  it  propagates.  However, 
if  a  portion  of  computation  proceeds  in  the  direction  opposite  to  that  of  wave 
propagation,  then  numerical  errors  may  grow  and  the  results  can  become 
meaningless,  which  is  likely  to  occur  when  the  path  of  propagation  is  long. 

The  dynamic  behavior  of  a  junction,  where  two  or  more  wave  guides  meet,  can 
also  be  characterized  by  a  scattering  matrix.  Figure  2  shows  a  junction, 
labeled  junction  A  and  possessing  of  several  ports,  each  of  which  represent  an 
interface  with  a  waveguide.  Each  interface  is  identified  by  a  number;  e.g. 
number  j  shown  in  the  figure.  Then,  according  to  the  convention  adopted 
earlier,  yj(A,j)  is  an  outgoing  wave  vector  and  Pj(j,A)  is  an  incoming  wave 
vector  with  respect  to  junction  A. 

A  governing  equation  for  junction  A  can  be  formulated  using  any  suitable 
method  of  structural  dynamics,  e.g.,  a  finite-element  procedure, 

w(l)  f(l) 


m  {  •  i  =  i  •  } 
•  • 
w(m)  f(m) 


8 


(18) 


where  [P]  is  an  impedance  matrix,  each  pair  of  w(j)  and  f(j)  are  displacement 
and  force  vectors,  respectively  at  interface  j  and  j=l,2,...,m.  The  state 
vector  at  the  jth  interface  can  be  transformed  into  wave  vectors  by 


^w(j) 

f(j) 


}  =  [D(j)]  { 


yj(A.j) 

yj(j.A) 


}  - 


Mj(A.j) 
Nj (A, j ) 


Mj(j.A)  Mj (A, j ) 

{  },  j=l 

Nj(j.A)J  Pj ( j ,A) 


,m,  (19) 


where  [Ou)3  is  constructed  from  the  eigenvectors  of  the  transfer  matrix 
associated  with  the  cell  belonging  to  the  waveguide  next  to  j.  Again  the  Mj  and 
Nj  submatrices  account  for  the  contributions  of  wave  vectors  toward 
displacements  and  forces,  respectively.  Substituting  (19)  into  (18)  yields 


Pi (A, 1 )  Pi (1 »A) 

{  •  1  =  CS]  {  .  }  (20) 

*  • 

Pm(A,m)  Pm(m,A) 


where  [S]  is  the  required  scattering  matrix  for  junction  A,  and  is  given  by 


[S]  =  (CP] 


Mi(A,l) 

Ni(A,l) 

0 

0 

0  ’. 

0  ’. 

Mm(A,m). 

Nm(A,m). 

Ml, A) 

Ml, A) 

.  0 

0 

• 

-  CP] 

• 

0 

0 

Nm(m,A). 

Mm(m,A). 

In  this  scattering  matrix,  a  total  of  m  square  submatrices  along  the  diagonal 
are  reflection  matrices,  reversing  incoming  waves  back  to  the  same  waveguides. 
The  remaining  m(m-l)  submatrices  are  transmission  matrices  redirecting  incoming 
waves  at  one  port  to  the  others.  In  the  special  case  of  m=2,  a  junction  reduces 
to  a  waveguide,  and  equation  (20)  reduces  to  equation  (8). 
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3.  A  WAVEGUIDE  SUBJECT  TO  AN  EXCITATION 


Figure  3  illustrates  a  waveguide  subject  to  an  excitation  at  j  which  is  the 
interface  between  cells  j-1  and  j.  The  excitation,  denoted  as  a  source  vector 
{ws,fs},  causes  a  discontinuity  in  the  wave  vector  at  the  interface,  given  by 


lij+(j-,j+)  Wj-(j-.j+) 

{  3  }  -  }  = 
Mj+(j+.j-)  Mj-(j+»j-) 


Mj(j-,j+)  Mj ( j- , j+) 
Nj(j-,j+)  Nj(j-,j+) 


-1  ws  X, 

{  }  -  L  J  (22) 

U  l* 


where  j-  and  j+  are  the  two  sides  of  interface  j,  and  Q(j  +  , j-)=[D(j)]_1[D(j-l)]. 
In  general,  only  one-half  of  the  elements  in  the  source  vector  are  non-zero. 
Equation  (22)  may  be  re-arranged  to  read 


yj-  (J+J-) 
Mj+  (J-J+) 


}  =  [S]  { 


(j-.  j+) 

(J+.  J-) 


}  + 


0 

I 


-  T ( j+,  j-)  ^  l2 

-  R(j+,  j-)  Zi 


(23) 


in  which  [S]  is  a  scattering  matrix  that  can  be  contructed  using  the  elements  of 
matrix  Q(j+,  j-). 

Now,  treating  the  two  sections  m  to  j-  and  j+  to  n,  as  two  separate 
waveguides,  equation  (23)  may  be  rewritten  as  follows: 


(J+.  J-) 

}  «  (I  -  [S] 

(j-,  j+) 


R(m.j-)  0 

0  R(j+,  n). 


([S] 

T(j-,m)  0 

fVm0n.j) 

l  }  + 

'  0 

-  T(j+,j-) 

0  T(n,j+)  . 

Mn(n,j) 

I 

-  R(j+J-)  . 

(24) 


Noticing  the  relationships, 

ymU.ra)  =  R(fn,j-)ym(m,j)  +  T(J-,m)yj_(J+,j-)  (25) 

Pn(j.n)  =  R(h, j+)pn(n»J)  +  T(j+,n)jij+(j-,j+)  (26) 

Mj+(j"»j+)  and  yj_(j+,j-)  can  be  eliminated  by  using  equation  (24),  arriving  at 
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R(m,n)  T(n,m) 
T(m,n)  R(n,m) 


^ym(m»n) j 
yn(n,m) 


yn(m,n) 


'T(j-.m)  O' 

* 

'R(J-.m)  0  ' 

-1 

'  0  -  T (j  +  ,j-) 

I  -  [S3 

.  I  -  R(j-.j+) 

0  T ( j+,n) 

0  R( j+,n) 

L 

L  }  (27) 

L 


in  which  the  identities  pm( j ,m)=ym(n,m) ,  etc.  have  been  used.  This  equation 
shows  how  the  relation  between  the  incoming  and  the  outgoing  waves  is  changed 
due  to  the  presence  of  external  excitation. 

Equation  (27)  is  written  in  a  form  indicative  of  the  adopted  computational 
scheme  in  which  the  outgoing  waves  are  treated  as  unknowns  and  solved  in  terms 
of  the  external  excitations  and  the  incoming  waves.  In  actuality,  the  incoming 
waves  are  also  unknowns,  and  two  additional  equations  are  needed  for  a  complete 
solution.  If  both  ends  of  a  waveguide  are  bounded,  then  the  required  two 
equations  can  be  obtained  in  the  form  of  equation  (13).  If  one  end  is  connected 
to  another  waveguide,  then  the  incoming  waves  for  one  waveguide  are  exactly  the 
outgoing  waves  for  the  connecting  waveguide.  Thus,  only  the  outgoing  waves  need 
be  treated  as  unknowns  and  numerical  stability  is  assured  in  the  computation. 


4.  STRUCTURAL  RESPONSE 

The  response  of  a  structure  at  a  certain  location  k  due  to  an  external 
excitation  at  another  location  j  can  be  formulated  most  simply  by  referring  to  a 
waveguide  which  contain  k  but  not  j.  This  arrangement  is  always  possible  as  can 
easily  be  explained  with  the  aid  of  figure  3.  Here  k  is  shown  to  be  between  m 
and  j.  Then  a  waveguide  extending  from  m  to  £,  where  £  is  between  k  and  j,  will 
meet  such  a  requirement.  Referring  to  this  waveguide,  we  obtain 


{ 


f  (k) 


Mk(m,k)  Mk(k,m) 

'CI-R(k.m)  R(M)]-1  0 

Nk(m,k)  Nk(k,m). 

0  [I-R(k,£)R(k,m)]'1 
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X 


yn(4,m) 


(28) 


T(m,k)  R(k,m)T(£,k) 

R(k,£)T'm,k)  T(£,k) 

It  is  of  interest  to  note  the  following  special  cases: 


4-1.  WAVEGUIDE  WITH  A  BOUNDED  END  AND  WITHOUT  EXTERNAL  EXCITATION 

In  this  case,  there  is  no  incoming  wave  at  the  bounded  end.  For  example, 
if  m  is  a  free  end,  then  ppm(Fm,£)=0.  The  displacement  response  at  k  can  be 
obtained  from  equation  (28)  as 

w(k)  =  [Mk(Fm,k)R(k,Fm)  +  Mk(k,Fm)][I-R(k,£)R(k,Fm)]-l  T(£,k)yn(£,k)  (29) 


4-2.  CLOSED  LOOP 

A  schematic  of  a  closed  loop  is  shown  in  figure  4.  An  external  excitation 
is  applied  at  a  cell-to-cell  interface,  labeled  1.  The  numbering  of  other 
interfaces  proceeds  clockwise,  assuming  that  an  observer  is  inside  the  closed 
loop.  The  excitation  gives  rise  to  a  discontinuity 


ws  w ( 1+)  w(l-) 


fs  ^(1+)  f(l-) 

For  simplicity,  let  the  two 
equation  (22) 


adjoining 


cells  at 


(30) 

be  identical.  Then  from 


ji1  +  (l+,l-)  ^Mi-(1+,1-J 

yl+(l-,l+)  Mi-(1-,1+) 


M,(l+,1-)  H1(l-,1+)1'1  ws^ 

Nxd+,1-)  N1(1-,1+)J  fs  L 


(31) 


To  render  equation  (28)  applicable,  the  closed  loop  is  treated  as  a  waveguide 
with  boundaries  at  1+  and  1-.  From  such  a  point  of  view,  the  waveguide  is 
without  an  external  excitation.  Clearly  then 


^1+(1-,1+)^  R(l+, 1-)  T(l-,1+) 

Mi-(1+,1-)  It(1+,1-)  R(l-,1+) 


^Mi+(1+»1-)^ 
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(32) 


Substituting  (32)  into  (31)  and  re-arranging. 


On  0 

'I 

-R(l-,1+) 

Pi-(1-,1+) 

.0  o22. 

R(l+,1-) 

I 

where 

o12-1  =  i  -  t ( i — , i+)  -  R(i+,i-)[i-T(i+,i-)]-iR(i-,i+) 
o22-‘  »  I  -  T ( 1+, 1-)  -  R(l-,l+)[I-T(l-,l+)]-lR(l+,l-) 


(33) 


(34a) 

(34b) 


Equation  (33)  describes  the  two  incoming  wave  vectors  at  the  two  ends  of  the 
waveguide  extending  clockwise  from  1+  to  1-.  Thus,  the  structural  response  can 
be  determined  using  equation  (28).  For  example,  the  displacement  response  at  k 
is  given  by 

w(k)  =  [Mk(l,k)  R(k , 1+)  +  Mk(k,l)][I-R(k,l-)R(k,l+)]-l  T(l-,k)y1_(l-,l+) 

+  [Mk ( 1 , k)  +  Mk(k,l)R(k,l-)][I-R(k,l+)R(k,l-)]-l  T(l+,k)yl+(l+,l-) 

k  *  1+,  2,  ...,  n,  1-  (35) 


5.  AN  ELIMINATION  PROCEDURE 

The  structural  response  discussed  in  Section  5  has  been  expressed  in  terms 
of  a  one-dimensional  waveguide.  To  obtain  the  response  for  more  complicated 
structures,  it  is  conceptional ly  convenient  to  convert  such  structures  to 
equivalent  one-dimensional  waveguides.  This  can  be  accomplished  by  use  of  an 
elimination  procedure  which  will  be  explained  below  for  various  cases  of 
increasing  complexity. 

5.1  A  WAVEGUIDE  WITH  AN  ADJOINING  JUNCTION  ON  ONE  END  AND  BOUNDED  ON  THE  OTHER 
Shown  in  figure  5  is  a  waveguide  which  interfaces  with  an  m-port  junction  A 
on  its  left  end  m  and  is  bounded  on  its  right  with  a  free  end  Fn.  Noticing  that 
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Vim(A,m)=pm(m,Fn)  and  nm(m,A)=ym(Fn,m) ,  we  obtain  from  substituting  (13)  into 

(20), 


Vx(A.l)  [I 


{  1  }  =  CS] 

PmMA.m-l) 


0 


0 


I 


Ml.A) 

{  •  } 
Vim-xtm-l.A) 


R(m,Fn) 


Mm(A.m) 


(36) 


which  is  valid  when  no  external  excitation  is  present  on  the  waveguide. 
Equation  (36)  can  be  simplified  to 


Vi(A,l)  Ml, A) 

•  • 

•  ru„  u12i  . 

I  •  }  -  I  _  I  (  •  I  (37) 

IM  u22j 

Ufi,-!  (a  ,m-l )  ym.1(m-l,A) 

o  ym(A,m) 

where  U^j  are  submatrices  of  the  product  of  the  two  square  matrices  on  the  left 
hand  side  of  equation  (36)  and  U22  =  U22  -  I.  Eliminating  wave  vector  pm(A,m) 
from  equation  (37),  we  obtain 

Vi(A.l)  Ml, A) 

(  }  =  [Sm-i]  {  .*  }  (38) 

•  • 

Pm-i(A.m-l)  Vm-i (m_l ,A) 

_-l 

in  which  Sm_i  =  Uu  -  U12  U22  U2l.  Equation  (38)  suggests  that  may  be 
considered  as  a  scattering  matrix  for  a  new  junction  with  m-1  ports. 

If  an  external  excitation  is  applied  on  a  bounded  waveguide,  then  equation 
(27)  must  be  employed  before  implementing  the  elimination  procedure,  in  which 
case  a  nonhomogeneous  term  due  to  the  excitation  will  appear  in  equation  (38). 
For  the  case  of  several  bounded  waveguides  intersecting  at  a  junction,  the  same 
elimination  procedure  can  be  applied  sequentially. 
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5.2  A  WAVEGUIDE  WITH  AN  ADJOINING  JUNCTION  ON  EACH  END 

As  shown  in  figure  6,  a  waveguide  is  connected  to  an  m-port  junction  A  at 
its  left  end  m,  and  connected  to  a  k  +  1-port  junction  B  at  its  right  end  n. 
External  excitations  are  applied  at  a  station  j  on  the  waveguide  between  m  and 
n.  To  proceed,  we  first  assemble  the  scattering  relationships  for  both 
junctions  as  follows 

Pi (A, 1 )  Pi (1 »A) 


fPm(A,m)} 

Pn(B,n) 


Sa  o 

(PmOM)! 
o  SBJ  Pn(R»B) 


(39) 


•  • 
pn+ic(B,n+k)  yn+k(n+k,B) 

where  and  SB  are  scattering  matrices  for  junctions  A  and  B,  respectively. 
Rearrange  equation  (39)  by  moving  wave  vectors  at  m  and  n  to  the  last  row  to 
obtain 

Pi(A.l)  Pi(l.A) 

•  • 

•  • 

{Pm(B,n+l)}  =  [S]  {pn+i(n+l,B)}  (40) 

•  • 

•  • 

•  • 

pm(A,m)  Pm(m>A) 

Pn(B.n)  Pn(n.B) 

Note  that  the  wave  vectors  at  the  two  ends  of  a  waveguide  can  be  related  in  the 
form  of  equation  (27)  and  that  PmCA.m^p^m.n) ,  pm(m,A)=pm(n,m) ,  pn(B,n)=pn(n,m) 
and  pn(n,B)=pn(m,n) .  Substitution  of  equation  (27)  into  equation  (40)  yields 
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VU(A,1) 


Vi(l.A) 


Mfji-i  (A,m-1) 

{yn+x (B.n+i) } 


ym-i(m_1»A) 

[Wn+i(n»®)  > 


pn+k(B,n+k) 

ym(A,m) 

Wu(B.n) 


R(m,n)  T(n,m) 
T(m,n)  R(n,m) 


yn+k(n+k,B) 

Vm(A.m) 

Un(B»n) 


+  [S]  {  .  } 


where 


rm  fT(j-,m)  0  1  r  fR(j-,j+)  T(J+,J-)1  fRCJ-.m)  0  1  1-1 

[  I  ■  I- 

rn  0  T(j+,m)J  U(j-,j+)  R(j+,j-)J  0  R(j+,n)J  . 


0  -T(j+J-)1  Y, 

x  {  }  (42) 

.  I  -R(j-J+  L 

Equation  (42)  implies  that  when  computing  structural  response  at  a  location 
outside  the  waveguide  m-n,  the  effects  of  the  external  excitations  at  j  can  be 
accounted  for  by  "equivalent"  excitations  at  stations  m  and  n. 


To  proceed  further,  denote 
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1(1) 

0 

0 

X(m-l) 

{X(n+1) }  =  [S]  {  .  }  (43) 

•  • 

rn 

X(n+k) 

X(m) 

X(n) 

Making  use  of  the  same  elimination  procedure  as  that  of  Section  6-1,  we  obtain 


Vi(A,l)  Ml, A)  X(l) 


ym-x(A,m-l)  __x  Pm-iOn-l.A)  Z(m-l)  _-l  J(m) 

{ym+i (B,n+1) }  *  [Ux i_U12U22U21]  {Pn+i (n+l »B) }  +  {Z(n+1)}  -  U12U22  {  } 

ZOO 


pn+k(B.n+k)  Mn+k(n+k,B)  X(n+k)  (44) 


in  which  U-jj  are  the  submatrices  of  the  product  of  the  two  square  matrices  on 
the  right  hand  side  of  equation  (41)  and  U22=U22-I.  It  is  of  interest  to  note 
that  although  the  wave  vectors  at  the  two  ends  of  the  waveguide  do  not  appear  in 
equation  (44),  their  effects  are  included  in  an  equivalent  new  junction  with 
ports  1,  2,  ...,  m-1  and  n+1,  ...,  n+k.  When  external  excitations  are  absent 
equation  (44)  reduces  to  a  relationship  among  the  scattering  waves  at  these 
ports. 


6.  COMPLEX  STRUCTURAL  NETWORKS 

Results  obtained  in  the  previous  sections  can  now  be  applied  to  complex 
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structural  networks.  To  illustrate,  consider  a  specific  structure  shown  in 
figure  7  in  which  waveguides  are  represented  by  solid  straight  lines  and  are 
labeled  with  circled  numbers,  and  junctions  are  denoted  by  capital  letters  A 
through  F.  An  excitation  is  applied  on  waveguide  (2)  and  the  response  is  to  be 
calculated  at  a  given  location  on  waveguide  (1).  A  number  without  circle 
indicates  the  interface  between  a  port  and  a  waveguide.  The  reason  for 
assigning  numbers  (1)  and  (2)  to  the  two  waveguides  where  the  output  and  the 
input  are  located  will  become  obvious  in  the  following  development. 

First,  we  assemble  the  wave  vectors  at  the  ports  of  all  junctions  and 
obtain  an  equation 


MA«#>  [Sa 

y#(B,#)  SB 

{  *  }  = 


y#(M) 

{  ‘  1 


MM) 


M#*F) 


(45) 


where  S/\,  SB,  ...,  and  Sp  are  scattering  matrices  for  junctions  A,  B,  ...,  and 
F,  respectively,  and  where  each  indicates  an  interface.  It  is  to  be 
understood  that,  when  representing  either  incoming  or  outgoing  waves  with 
respect  to  a  multi-port  junction,  the  symbol  It  ranges  over  all  the  ports.  For 
instance,  referring  to  figure  7,  p#(A,#)  represents  all  the  outgoing  waves  from 
junction  A;  namely,  it  is  an  abbreviation  of  the  vectors  (A, 1) } ,  {p,(A,8)}, 
and  {iin(A,U)}. 

Rearrange  the  order  of  wave  vectors  In  equation  (45)  in  the  numerical  order 
of  the  interfaces 


MM)  Ml,*) 

{  .  !  ■  [S]  {  }  (46) 

y2o(*,20)  y2  o (20 ,*) 

where  each  signifies  a  particular  junction  from  A  to  F.  The  new  scattering 
matrix  [S]  In  equation  (46)  is  of  course  not  a  band  matrix.  Now,  for  the 
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waveguides  (3)  to  (11),  the  scattering  relationship  may  be  written  as 


=  LSkJy#(*,#)  k  =  3,  11  (47) 

in  which  each  represents  one  particular  end  of  the  kth  waveguide  and 
indicates  the  junction  to  which  the  kth  waveguide  is  connected.  Moreover,  as 
shown  in  figure  7,  (3)  and  (11)  are  bounded  on  one  end.  Thus,  S2  and  S1X  are 
reflection  matrices  of  order  pxp  instead  of  2px2p.  Substituting  (47)  into  (46), 
we  obtain 


Mi(*,l) 

• 

10 

Mi (1 ,*) 

• 

M* (*,4) 

}  =  [S] 

Ms(*,5) 

s. 

M*(4>) 

(Ms(*,5)} 

• 

• 

• 

0 

Mj  o (*,20) 

Bn 

Mjo(*,20) 

(48) 


The  same  elimination  procedure  employed  before  can  be  applied  to  eliminate 
ys  through  y20  in  equation  (48),  resulting  in 


Pi(A,l) 

V*(B,2) 

y,(D,3) 

MF,4) 


}  -  C5]  f 


Ml, A) 

M2,B)^ 


V»(3,D) 

M*(4,F) 


(49) 


Since  the  excitation  is  applied  on  waveguide  (2),  the  scattering  matrix  for 
waveguide  (2)  can  be  used  to  reduce  equation  (49)  to 


Mi(A.l) 
MB, 2) 


,  _-l  Mi ( 1 » A) 

}  =  [Ml  l-Uj  2  1^2  2  1^2  1  U  { 

Mi (2,B) 


i . 

1(2) 


_ ~ 1  .ZO), 

MM  {  } 

1(4) 


(50) 


according  to  equation  (44).  Equation  (50)  Indicates  an  input-output 
relationship  for  a  one-dimentional  waveguide  which  is  equivalent  to  the  network 
containing  all  the  intersections  and  waveguides  (2)  through  (11).  Noticing  that 
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V»I  (1  .A)  Vi(A,l) 

{  }  -  [St]  {  } 

y*(2,B)  MB, 2) 


(51) 


the  incoming  waves  at  the  two  ends  of  waveguide  (1)  are 

Hi(A,l)  ,  Id),  _lfE(3), 

{  }  =  [I-tUn-U^Uliu.JCS,]]  1  li  }  -  U12UlJ  {  (52) 

lit(B,2)  1(2)  1(4) 

Finally,  the  response  at  an  arbitrary  location  on  waveguide  (1)  can  be 

calculated  from  equation  (28). 

Notice  that  a  matrix  inversion  is  required  when  obtaining  equation  (49) 
from  equation  (48).  This  inversion  can  be  expensive  when  the  number  of 

waveguides  is  large.  Alternatively,  instead  of  assembling  all  the  junctions  in 
the  beginning  as  shown  in  equation  (45),  one  can  carry  out  the  elimination 

procedure  one  waveguide  at  a  time,  starting  with  one  or  two  junctions.  More 
junctions  are  added  successively,  until  the  entire  process  is  complete.  Such  a 
step-by-step  elimination  scheme  involves  inversions  of  smaller  matrices,  and 
requires  much  less  computational  efforts.  For  the  specific  example  shown  in 

figure  7,  we  can  begin  with  junction  F,  and  eliminate  first  waveguide  (11).  Add 
successively  junctions  E,  D,  and  C,  and  eliminate  waveguides  (10),  (9),  through 
(3)  one  by  one. 

7.  A  NUMERICAL  EXAMPLE 

Numerical  calculation  has  been  carried  out  for  a  hypothetical  space 
structure  shown  in  figure  8.  The  global  form  of  the  structure  has  two  closed 
loops,  six  junctions  and  two  free-ended  arrays.  The  array  in  the  N-S  direction 
shared  by  the  two  closed  loops  consists  of  32  cells.  Two  additional  N-S  arrays 
with  32  cells  each  and  four  E-W  arrays  with  41  cells  each  complete  the  two 
closed  loops.  The  free-ended  array  to  the  north  are  composed  of  40  cells  and 
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the  free-ended  array  to  the  south  composed  of  80  cells.  The  two  cells,  each  of 
which  is  shared  by  four  intersecting  arrays  are  treated  as  four-port  junctions. 
The  four  corner  cells  are  treated  as  two-port  junctions. 

The  construction  of  a  typical  cell  unit  is  depicted  in  figure  9.  The  two 
diagonal  bars  on  each  face  are  assumed  to  be  unconnected  at  the  intersection. 
All  bars  are  of  a  cross-section  area  of  5.02  x  10'3  m2 ,  and  are  pinned  at  the 
corners  of  the  cell;  thus,  bending  moments  cannot  be  developed  in  the  bars.  In 
this  case,  the  structure  cannot  resist  the  rotational  degrees  of  freedom  at  the 
corners,  and  only  three  translational  degrees  of  freedom  per  corner  need  be 
retained  in  the  analysis.  Still,  for  this  model  structure,  the  total  number  of 
degrees  of  freedom  is  over  5000. 

The  material  properties  used  in  the  calculation  include  Young's  modulus  E  = 
69.6  x  109  (1+iy  sgn  w)  N/m2,  where  Y  =  0.01,  and  the  mass  density  p  =  3.77  x 
103  kg/m3. 

Referring  to  figure  8,  an  excitation  in  the  E-W  direction  is  applied  at 
station  j  which  is  at  the  middle  of  an  N-S  array.  The  corresponding  force 
responses  are  calculated  at  two  stations:  station  h  and  station  k.  Station  k 
is  located  at  the  interface  between  the  16th  and  the  17th  cells  counting  from 
the  western  end  of  an  E-W  array  as  shown,  while  station  h  is  at  the  interface 
between  the  10th  and  the  11th  cells  counting  from  the  southern  free  end  of  an  N- 
s  array. 

Shown  in  figure  10  is  the  computed  magnitude  of  the  frequency  response 
function  HjJw)  for  the  vertical  shear  force  response  at  k.  Significant  response 
Is  found  for  frequencies  as  high  as  550  rad/s,  apparently  involving  many  modes. 
Thus,  the  usual  assumption  of  dominant  lower  modes  in  a  conventional  modal 
analysis  can  lead  to  significant  errors  in  this  case.  Figure  11  shows  the 
corresponding  Impulse  response  in  the  time  domain.  No  response  is  felt  at  k 
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until  about  0.0388  seconds  when  fastest  traveling  waves  arrive.  These  waves  are 
associated  with  higher  frequencies,  and  they  travel  through  the  shortest  path 
containing  77  cells.  The  time  delay  of  0.0388  suggests  that  the  speed  of  the 
first  arriving  wave  is  about  2000  cells  per  second.  After  reaching  a  maximum  at 
t  =  0.048  sec,  the  response  gradually  dies  down  until  about  0.666  sec,  when  the 
second  wave  group  arrives  also  through  the  shortest  path.  The  response  is 
reinforced  further  by  the  higher  velocity  waves  propagating  through  the  second 
shortest  path  containing  144  cells.  The  very  complicated  time-history  after 
0.095  sec  indicates  the  superposition  of  waves  traveling  through  different  paths 
at  different  velocities  and  scattering  at  the  junctions. 

The  frequency  response  and  the  corresponding  impulse  response  of  the 
induced  horizontal  shear  force  at  h  are  shown  in  figures  12  and  13, 
respectively.  It  is  interesting  to  note  that  in  contrast  to  the  case  of 
location  k,  lower  modes  with  frequencies  ranging  from  0  to  200  rad/s,  now 
dominate  the  response.  The  shortest  path  from  j  to  h  contains  127  cells;  thus 
It  takes  0.064  sec  for  the  fastest  traveling  waves  to  arrive  to  the  site,  as 
indicated  in  the  impulse  response  plot. 

8.  CONCLUDING  REMARKS 

In  the  above  analysis,  attention  is  focused  on  structural  response  at  the 
interfaces  of  truss  units,  bypassing  detailed  motions  in  the  individual  bars 
within  each  unit.  This  is  in  line  with  the  usual  treatment  of  wave  propagation 
in  periodic  structures  [4].  If  responses  of  certain  Individual  bars  are  of 
interest,  then  these  bars  and  perhaps  some  otner  surrounding  bars  must  be 
treated  also  as  waveguides.  In  such  a  case,  the  theory  developed  herein  is 
still  applicable,  except  for  the  additional  expense  in  computer  time.  The  need 
to  investigate  individual  structural  members  appears  quite  realistic  for  civil 
engineering  structures. 
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Fig.l  A  waveguide  with  incoming  and  outgoing  waves  at  the  ends. 
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Fig. 2  A  m-port  junction. 
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Fig. 10  Amplitude  of  frequency  response  of  a  shear  force  at  station 
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Fig. 13  Impulse  response  of  a  shear  force  at  station 


LOCALIZATION  OF  WAVE  PROPAGATION  IN  DISORDERED  PERIODIC  STRUCTURES 

Due  to  material,  geometric  and  manufacturing  variabilities,  a  structure 
designed  to  be  periodic  can  never  be  exactly  periodic.  The  departure  from  ideal 
periodicity  is  known  as  disorder,  which  must  be  taken  into  account  in  a 
realistic  analysis. 

The  ability  to  transmit  disturbances  indefinitely  within  the  wave-passage 
bands  is  reduced  when  the  ideal  periodicity  in  a  periodic  structure  is 
disrupted.  This  is  known  as  the  localization  effect,  first  predicted  by 
Anderson  in  a  celebrated  paper11  concerning  the  transport  of  electrons  in  an 
atomic  lattice,  which  is  essentially  a  disordered  three-dimensional  periodic 
system.  Anderson's  work  has  found  important  applications  in  solid  state 
physics.  Since  then,  numerous  papers  on  localization  have  appeared  in  the 
physics  literature12. 

Obviously,  localization  also  occurs  in  disordered  periodic  structural 
systems,  but  Its  theoretical  investigation  is  more  difficult  than  that  of  a  one- 
dimensional  atomic  lattice  since  governing  equations  for  structural  systems  are 
generally  more  complex.  Using  the  models  of  a  chain  of  coupled  pendula  and  a 
vibrating  string  constrained  by  attached  masses  and  springs,  Hodges13  provided 
an  excellent  explanation  for  the  confinement  of  vibration  and  for  localized 
normal  modes  as  a  result  of  structural  irregularity.  Subsequently,  he  and 
Woodhouse14  calculated  the  so-called  localization  factor  for  one  of  the 
systems,  characterizing  the  average  exponential  decay  of  vibration  amplitude 
with  respect  to  the  distance  from  the  disturbance.  Two  perturbation  schemes 
were  devised  to  treat  the  cases  of  weak  and  strong  localizations,  depending  on 
the  relative  measures  of  disorder  and  internal  (namely,  unit-to-unit)  coupling. 

Pierre1^  also  calculated  the  localization  factor  for  both  weak  and  strong 
localization,  using  both  the  wave  propagation  and  normal  mode  formulations  and 
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obtained  essentially  same  results.  He  devised  two  perturbation  schemes  to  treat 
the  two  cases  when  the  ratio  of  the  disorder  measure  to  the  internal  coupling 
measure  is  either  large  or  small.  As  expected,  his  results  were  not  very 
accurate  when  this  ratio  is  of  order  one. 

The  first  analysis  of  a  generic  disordered  periodic  system  was  given  by 
Kissel16,  making  use  of  the  concept  of  wave  transmission  and  reflection,  and  a 
limit  theorem  on  the  product  of  random  matrices  due  to  Furstenberg1?.  For  a 
wave  motion  propagating  at  a  wave-passage  frequency,  no  reflection  would  occur 
at  the  unit-to-unit  interface  if  the  system  were  exactly  periodic.  When  two 
adjoining  units  are  not  identical,  total  transmission  is  no  longer  possible,  and 
some  reflection  will  take  place  at  the  interface.  However,  Kissel's  analysis 
requires  that  the  reflection  be  small;  thus,  his  results  are  valid  only  for  the 
case  of  weak  localization. 

Unlike  previous  investigators,  Kissel  has  also  explored  the  case  of  multi¬ 
coupling;  namely,  the  internal  coupling  which  permits  the  transmission  of 
several  types  of  waves.  An  outline  was  provided  of  an  approach  based  on  the 
multiplicative  ergodic  theorem  of  Oseledets1?.  Unfoi tunateiy,  the  theorem  only 
predicts  the  existence  of  a  localization  factor  (analogous  to  the  Lyapunov 
exponent  in  the  time  domain,  of  interest  to  Oseledets),  without  providing  the 
necessary  clues  as  how  it  can  be  calculated. 

The  ability  to  treat  a  generic  system  is  clearly  very  desirable,  since  the 
analysis  is  applicable  to  a  general  class  of  problems,  not  restricted  to  a 
specific  type  of  governing  equations.  Furthermore,  it  appears  easier  to  extend 
such  an  analysis  to  the  case  of  multi-coupling  which  is  likely  to  occur  in  many 
practical  structural  systems.  Thus,  it  is  also  the  objective  of  the  present 
research  to  calculate  the  localization  factor  for  a  generic  disordered  periodic 
structure.  A  new  perturbation  scheme  is  developed  here  which  permits  successive 
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improvement  of  the  accuracy,  making  it  applicable  to  either  weak,  moderate  or 
strong  localization.  The  additional  accuracy  is  obtained  by  taking  into  account 
reflections  from  additional  nearby  cells.  Account  is  also  taken  of  structural 
damping;  thus,  the  results  are  more  useful  to  practicing  engineers.  The  key  to 
the  successful  development  of  a  simple  theory  under  rather  general  settings  is 
the  substitution  of  ensemble  averages  for  sequential  averages  of  certain 
statistical  properties  of  the  structure  on  the  basis  of  spatial  ergodicity.  The 
simplicity  of  the  new  method  is  illustrated  by  applying  it  to  a  multiply 
supported  Euler-Bernoul  1  i  beam  with  an  additional  torsional  spring  at  each 
support.  The  accuracy  of  the  results  is  substantiated  by  Monte-Carlo 
simulations.  Only  mono-coupling  at  the  unit-to-unit  interfaces  is  considered  in 
detail.  The  case  of  multicoupling  will  be  considered  in  the  future. 


Wave  Transmission  and  Reflection  Matrices 
Consider  two  typical  cell-units  in  a  disordered  periodic  structure,  denoted 
as  cells  n  and  n+1  and  shown  in  Fig.  la.  The  relation  between  the  state  vector 
at  the  left  end  of  cell-unit  n  and  that  at  the  left  end  of  cell-unit  n+1  may  be 
written  as  follows 


w(n+l£)  w(n£) 

f  J  -  CT ( n) 3  {  } 

f(n+l i)  f(n£) 


(1) 


where  each  w  is  a  p-dimensional  vector  of  generalized  displacements,  each  f  is  a 
p-dimensional  vector  of  generalized  forces,  and  [T(n)]  is  a  2px2p  transfer 
matrix  associated  with  the  nth  cell.  It  is  implied  in  Eq.(l)  that  the  motion  is 
time-wise  sinusoidal,  and  the  displacements  and  forces  are  interpreted  as  their 
complex-valued  amplitudes.  Because  of  disorder,  [T( n ) ]  is  generally  different 
from  the  ideal  transfer  matrix  [T]  associated  with  the  ideal  design  condition  of 
exact  periodicity.  Nevertheless,  it  is  always  possible  to  write 
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CT(n)]  =  [P(n)][T]  (2) 
which  means  that  the  effect  of  disorder  in  the  nth  cell  can  be  lumped  at  its 
right  end  and  represented  by  a  "point"  transfer  matrix  [P(n)].  In  other  words, 
the  original  system  shown  in  Fig.  la  can  be  replaced  by  an  equivalent  system 
shown  in  Fig. lb.  In  typical  transfer  matrix  analysis19,  a  point  transfer  matrix 
is  used  to  relate  state  vectors  on  the  two  sides  of  a  structural  discontinuity 
which  may  arise  from  a  concentrated  mass,  spring,  damper,  etc.  However,  [P(n)] 
in  Eq. (2)  is  merely  a  mathematical  device,  and  it  may  even  correspond  to,  e.g., 
a  negative  mass. 

Refer  to  Fig.  lb  and  write 


w(nr)  w(n£) 

J  =  [T]  {  } 

f ( n  r )  f(n£) 

(3-a) 

w(n+l£)  w(nr) , 

1  -  CP(n)]  {  } 

f(n+l£)  f(nr) 

(3-b) 

Each  state  vector  in  Eq.(3-a)  can  be  transformed  into  a  wave  vector  using 


w(n4)  Mr(n£) 

{  1  -  IM  .  ! 

f(" l) 


(4-a) 


w(nr)  Mr(nr) 

{  1  ■  CD]{  o  } 

f(nr)  y*(nr) 


(4-b) 


where  the  columns  in  the  transformation  matrix  [D]  are  the  eigenvectors  of  [T]. 
These  eigenvectors  must  be  arranged  in  a  particular  order,  which  will  become 
apparent  in  the  subsequent  development.  We  then  obtain  a  wave  transfer 
relationship 


Vr(nr)  ,  ryr(na) 

{  )  -  CD]-1[T][D]{  }  = 


^Vr(n£)^ 


(5) 
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y£(nr)  y£(na)  I  .  I  y£(n£) 

[  0  ^  2  pi 

in  which  A j ,  j  =  l,  2p  are  the  eigenvalues  of  the  ideal  transfer  matrix  [T]. 
Each  Aj  represents  the  ratio  of  complex-valued  amplitudes  of  one  type  of  wave  at 
the  two  ends  of  an  ideal  cell-unit.  Since  damping  generally  exists  in  a 
physical  structure,  the  magnitude  of  a  wave  decreases  when  propagating  across  a 
cell.  Thus,  each  eigenvalue  accounts  for  the  reduction  of  the  wave  magnitude, 
as  well  as  the  change  of  phase  from  one  side  of  an  ideal  cell-unit  to  another. 
It  is  known  that  eigenvalues  of  a  transfer  matrix  are  reciprocal  pairs^.  For 
every  A j ,  there  exists  another  A|<  =  Aj-1.  If  |  Aj  |  <  1,  then  |  A|<  |  =  |  Aj  |  _1  >  1. 
They  correspond  to  two  identical  types  of  waves,  except  that  one  is  right-going 
and  the  other  is  left-going.  The  limiting  case  |  Aj  |  =  1  can  only  occur  in  an 
undamped  structure,  and  only  at  a  frequency  in  one  of  the  wave-passage  frequency 
bands. 

We  arrange  the  order  of  eigenvectors  in  [D]  such  that  the  corresponding 
eigenvalues  are  in  an  ascending  order;  thus  |  Aj  |  <  1,  if  j  =  1,  2,  ...,  p  and 
|  Aj  |  >  1,  if  j  =  p+1,  p+2,  ...,  2p.  Eq.  (5)  may  then  be  re-written  as  follows 


^yr(nr)  A 

y£(nr)  [o 


0  yr(n£) 

-1  l  „  J 

A  y£(n£) 


where 


A 


Ax  0 

0  Ap 


(6) 


(7) 


Or,  alternatively. 


^£("i) 

yr("r) 


}  -  [S]{ 


yr("a)j 

y£(nr) 


(8) 
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with 


The  meanings  of  the  superscripts  r  and  l  associated  with  the  wave  vectors  y  are 
now  clear;  they  indicate  the  right-going  and  the  left-going  wave  groups, 
respectively,  transmitted  from  one  end  of  an  ideal  cell  to  the  other,  as  shown 
schematically  in  Fig.  2.  In  fact,  Eq.(8)  has  been  written  to  suggest  that 
yr(n£)  and  y£(nr)  are  the  input  waves  and  y*(n£)  and  yr(nr)  are  output  waves,  as 
viewed  by  the  structural  unit  which  is  characterized  by  [S]. 

Similar  transformation  of  the  state  vectors  in  Eq.(3-b)  to  wave  vectors 
results  in 


^yr(n+l£) 

y*(n+l£) 


}  =  [Q(n)]{ 


yr(nr) 

y*(nr) 


(10) 
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In  Eqs.  ( 13a- 13d )  Sij  and  qij  are  pxp  submatrices  of  [S(n)]  and  [Q(n)], 
respectively,  and  the  argument  (n)  is  omitted  in  denoting  the  submatrices  for 
simplicity. 

It  is  physically  meaningful  to  re-name  the  submatrices  of  [$(n)]  as 
rr(n)=s11 ,  t£(n)=s12,  tr(n)=s21,  and  ri(n)=s22,  and  re-write  Eq. (12)  as 

y£(nr)  [rr(n)  t£(n) 

{  }  = 

yr(n+l4)  _tr(n)  r£(n) 

The  submatrices  t  and  r  are  called  the  transmission  and  reflection  matrices, 
respectively.  The  superscript  for  each  of  these  submatrices  signifies  the 
direction  of  an  incoming  wave  group  to  be  transmitted  or  reflected.  For 
example,  tr(n)  is  a  transmission  matrix  for  an  incoming  wave  group  traveling  in 
the  right  direction. 

Matrix  [S(n)],  which  contains  the  transmission  and  reflection  submatrices, 
is  known  as  a  scattering  matrix.  This  particular  scattering  matrix 
characterizes  the  effects  of  disorder  in  cell  n  which  are  lumped  at  the  right 
end  of  cell  n,  as  shown  schematically  in  Fig. 2.  In  wave  mechanics,  the  concept 
of  a  scattering  matrix  is  quite  general;  it  characterizes  the  behavior  of  any 
identifiable  structural  element.  In  fact,  matrix  [S]  in  Eq . (8)  is  also  a 
scattering  matrix,  in  which  the  transmission  and  reflection  matrices  are 
tr =t£=A  and  rr=r£=0.  These  submatrices  of  [S]  are  also  related  to  those  of  a 
matrix  [Q]  according  to  Eqs. (13a-d) ,  although  in  this  special  case 


It  is  interesting  to  remark  that  all  scattering  matrices  are  symmetric,  a 
consequence  of  the  reciprocity  theorem  for  linear  acoustic  systems. 

For  a  chain  of  cell  units,  1  through  N,  we  may  write 


jMr(nr)  j 

y2(n+la) 
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(16) 


yr(N+U) 

j|4(N+lA) 


}  =  [Q(N,1)]{ 


Mr(U) 

y4(i*) 


} 


where 

CQ(N,1)]=[Q(N)][Q]CQ(N-1)][Q]...[Q(1)][Q]  (17) 

and  where  [Q(n)]  and  [Q]  are  given  in  Eqs.  (11)  and  (15),  respectively.  Eq.(16) 
can  be  re-cast  in  the  form  of 


Va(1a) 

yr(N+l£) 


}  =  [S(N  *  1 ) ]{ 


yrd£)  } 

y4(N+l£) 


in  which 

[S(N, 1)] 


rr(N,l)  t*(N,l) 
tr(N, 1)  ra(N,l) 


(18) 


(19) 


Elements  of  matrix  [S(N,1)]  are  related  to  those  of  matrix  [Q(N,1)]  according  to 
Eqs.  (13a-d).  From  the  second  row  of  Eq.(18),  we  obtain 

yr(N+l£)  =  tr(N,l)yf(i£)  +  r*(N,l)y*(N+l£)  (20) 


The  Case  of  Mono-Coupling 

We  shall  now  consider  the  case  of  mono-coupling  in  more  detail.  In  this 
case  the  wave  transfer  matrices  Q  and  the  scattering  matrices  S  become  2x2,  and 
the  transmission  and  reflection  matrices  reduce  to  the  transmission  and 
reflection  coefficients.  Moreover,  we  shall  assume  for  the  time  being  that 
cells  N+l ,  N+2,  ...,  •  are  without  disorder;  thus  yfi(N+l£)=0.  Then  Eq.(20)  is 
simplified  to 


yr(N+l*)  =  tr(N,l)yr(li) 

The  localization  factor  is  defined  as 


(21) 
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(22) 


Y  =  -  Aim 
N— 

y£(N+ia) 


Hr(N+lg) 

yrdi) 


Upon  using  Eq.(21) 

Y  =  -  Aim  -  An  |  tr(N, 1)  |  (23) 

N-*»  n 

It  implies  that,  if  the  above  limit  exists,  the  transmitted  wave  has  an  average 
exponential  decaying  rate  of  y  per  cell-unit.  This  simple  result  can  also  be 
derived  from  an  application  of  Furstenberg' s  theorem  on  products  of  random 
matrices17. 

In  passing,  we  remark  that,  for  the  mono-coupling  case,  tr(N,l)  =  tA(N,l) 
as  can  be  proved  from  Eqs.(13a-d)  and  from  the  fact  that  |  Q(N,1)  |=  1  since 
[Q(N,1)]  is  a  transfer  matrix3. 

Several  specific  cases  will  now  be  considered  in  the  order  of  increasing 
complexity: 

(1)  Small  reflections  —  If  all  the  reflection  coefficients  due  to 
disorder  are  small,  such  that  the  effect  of  multiple  reflections  is  negligible, 
then  tr(N, 1)  can  be  approximated  by 

tr(N,l)  =  tr(l)  tr(2)  ...  tr(N)AN  (24) 

From  Eq.  (13b),  tr(n)=t£(n)=t(n)=q22_1 (n) .  Whence 

i  N 

Y  =  -  An  |  A  |  +  Aim  ±  l  An  |  q22(n)  |  (25) 

N-*»  N  n=l 

For  disordered  cells,  each  uncertain  parameter  may  be  represented  as  a  mean  (or 
nominal  value)  plus  a  random  variable.  Then  q22(n)  are  functions  of  such  random 
variables  which  may  be  denoted  by 


q22(n)  =  q22[X(n)] 


(26) 
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where  X(n)  =  {Xx (n) ,X2 (n) . Xk(n) }  is  a  vector  of  random  variables  with  zero 

means.  It  is  reasonable  to  assume  that  Xj(n)  for  different  n  are  independent 
and  identically  distributed,  and  that  they  form  an  ergodic  sequence.  Then  the 
sequential  averages  in  Eq. (25)  may  be  replaced  by  ensemble  averages;  namely, 

Y  =  -  an  |  X  |  +/  An  |  q22(x)  |  p(x)  dx  (27) 

where  p(x)  is  the  probability  density  of  X(n)  and  the  integrations  are  carried 
out  in  the  domain  of  the  random  vector  X(n).  Eqs.(25)  and  (27)  are  equivalent 

with  probability  one.  If  X(n)  is  a  k-vector,  the  integral  in  Eq.  (27)  is 
actually  a  k-fold  integral.  Given  the  forms  of  q22[X(n)]  and  the  probability 
distribution  of  X(n),  the  integral  can  be  evaluated  at  least  numerically.  In 
simple  cases,  it  can  often  be  carried  out  in  closed  form,  as  shown  later  in  an 
example. 

If  the  system  disorder  is  also  small,  then  Eq. (27)  may  be  approximated  by 
2 

°1  92 

Y  *  -  *n  |  A  |  +  7”  ^  r~  an  I  922(x)  |  }  x  =  0  (28) 

2  a*i 

where  o^2  are  the  variances  of  Xi(n).  If,  in  addition,  the  structure  is 
undamped  and  y  is  evaluated  at  a  wave-passage  frequency,  then  |  A  |  = 1 ,  and  Eq. 
(28)  reduces  to 

Oi2  32 

Y  =  ~  (  ~^n  I  Q 2 2  ( x )  |  } x= 0  (29) 

2  ax2 

Eq. (29)  has  been  obtained  previously  by  Kissel16. 

It  should  be  emphasized  that  Eq. (27)  is  valid  only  for  small  reflection 
coefficients  which  imply  that  localization  is  weak.  Eq.(28)  requires  both  small 
reflection  coefficients  and  small  disorder,  and  Eq. (29)  requires  additionally  a 
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zero  damping  and  a  wave-passage  frequency. 

(2)  Moderate  reflections  --  If  the  reflection  coefficients  due  to  disorder 
are  not  sufficiently  small,  the  approximation  Eq.(24)  may  not  be  adequate.  For 
improved  accuracy,  certain  multiple  reflections  will  be  considered.  As  shown  in 
Fig.  3,  the  right-going  wave  pr(nr)  gives  rise  to  a  right-going  transmitted  wave 
t(n)pr(nr)  and  a  left-going  reflected  wave  r(n)yr(nr).  The  reflected  wave 
travels  through  cell  n  and  is  itself  partially  reflected  back  at  the  interface 
between  cells  n  and  n-1.  This  back-reflected  wave  travels  also  through  cell  n, 
whereupon  it  splits  again  into  reflected  and  the  transmitted  parts.  The  latter 
is  clearly  t(n)A2r4(n-l)rr(n)yr(nr) ,  and  the  former  travels  back  through  cell  n, 
etc.  Summing  up  all  the  right-going  waves  which  pass  through  the  interface 
between  cells  n  and  n+1,  we  obtain 

pr(n+l4)=t(n){l+A2  ra(n-l)rr(n)+[A2  ra(n-l) rr(n)]2+. . . }yr(nr) 


1 


linl _ 

r£(n-l)  rr(n) 


yr(nr) 


(30) 


Thus 


N 

t r  (  n  ,  i )  =  r  i 


n«l 


t(n)A _ 

1  -  A2  r*(n-l)  rr(n) 


(31) 


It  follows  from  Eqs.(23)  and  (13a-d),  and  the  ergodicity  assumption. 


Y  =  -in  |  A  |  +  /  £n  |  q22(x)  |  p(x)dx 


+  //  in  |  1  + 


^  Qi a (x2 )q* i (Xj ) 
q2 2  (x2 )  q2 2 (Xx ) 


I  P(Xj)p(x2)  dx2dx2 


(32) 


The  third  term  Is  seen  to  provide  added  accuracy  over  Eq.(27).  Eq.(32)  is  valid 

for  the  case  of  moderate  localization. 
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(3)  Strong  reflection  —  Eq. (32)  may  still  be  Inadequate  if  the  reflection 
coefficients  are  relatively  large  (although  each  cannot  have  a  magnitude  greater 
than  one).  Higher  accuracies,  however,  can  always  be  obtained  by  including  more 
reflection  terms  in  the  formulation.  It  can  be  shown  that  if  the  effect  of 
reflection  at  the  interface  between  cells  n-1  and  n-2  is  also  taken  into 
consideration,  then 
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An  Example 


For  illustration,  the  above  theory  is  applied  to  an  Euler-Bernoull i  beam  on 
evenly  spaced  hinge  supports  and  with  an  additional  torsional  spring  at  each 
support,  as  shown  in  Fig.  4.  In  order  to  focus  our  attention  on  certain  key 
issues  in  the  theory  without  being  obscured  by  unnecessary  complexities,  it  is 
assumed  that  the  beam  is  undamped  and  that  only  the  torsional  spring  stiffnesses 
are  random  and  are  described  by 

kn  =  k0[l+  X(n)]  n  =  1,  2,  ...  (37) 

where  k0  is  the  average  kn  and  X(n)  are  random  variables  with  zero  means.  The 
other  physical  parameters  are  assumed  to  be  deterministic,  including  the 
distances  between  neighboring  supports  £,  the  bending  rigidity  of  the  beam  El, 
and  the  mass  of  the  beam  per  unit  length  m.  The  theory  is,  of  course, 
applicable  when  several  or  all  the  parameters  are  random. 

Several  choices  can  be  made  for  a  typical  cell-unit  in  the  analysis.  The 
one  selected  for  the  following  discussion  is  a  typical  beam  element  between  two 
neighboring  supports  plus  the  entire  right  spring.  The  entire  left  spring  is 
treated  as  belonging  to  the  preceeding  cell.  Then  the  transfer  matrix  for  the 
nth  cell  is  given  by3. 

a_ 
v 

+  |  [1+X(n)]  5  +  +  X(n)3 

2 

where  S  *  — ,  v  =  w  =  frequency 

£k  o  El 

coshv  cosv  -  1 

a  =  -  (39a) 

sinhv  -  sinv 


0 

[T ( n ) ]  = 
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and 


sinhv  t osv  -  coshv  sinv 

6  - -  (39b) 

sinh  v  -  sin  v 

The  non-dimensional  S  is  called  the  internal  (cell -to-cel 1 )  coupling  parameter, 
which  is  a  measure  of  relative  resistance  against  a  rotation  at  a  support 
between  the  beam  element  and  an  average  torsional  spring.  If  k0  is  very  large 
such  that  S  is  very  small,  then  there  is  little  ceil-to-cell  coupling,  and 
strong  localization  is  expected  to  occur19. 

The  transfer  matrix  [T]  for  the  ideal  cell-unit  without  disorder  and  the 
point  transfer  matrix  [P(n)]  lumping  the  effect  of  disorder  are  obtained  as 
follows 


[T] 


a 

v 


-  y(i  -j!)  + 1 


(40-a) 


1  0 


[P(n)] 


xin) 

S 


1 


(40-b) 


The  eigenvalues  of  transfer  matrix  [T]  are  a  reciprocal  pair  which  may  be 
denoted  conveniently  by  e±i9.  It  can  be  shown  that  0  can  be  found  from3 


cos0  =  0  +  — 

2Sv 

and  the  wave-passage  frequency  bands  are  determined  by  the  inequality 


(41) 


-15  0  + 


a 

26v 


5  1 


(42) 
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The  wave  transfer  matrix  [Q(n)]  is  obtained  from 


[Q(n)J 


l-iAX(n) 

iAX(n) 


-  iAX(n) 

1  +  i AX ( n ) 


(43) 


where  A  =  a/(26vsin8). 

The  following  physical  properties  have  been  used  in  the  numerical 
calculations:  l  =  0.1651m,  m  =  1.8043  kg/m,  and  El  =  0.3143  N-m2.  In  addition, 
X (n)  are  assumed  to  be  identically  and  uniformly  distributed  between  -  /3o  and 
/3o  where  o  is  the  standard  deviation  of  X(n). 

Substituting  Eq . (43 )  into  Eq. (27) ,  we  obtain  an  approximate  localization 
factor  for  the  case  of  wave-passage  frequency  and  weak  localization. 


Y  =  -£n(l+3A2c2)-l  +  — -  tan"1  (/3Ao)  (44) 

2  /3AO 

where  o2  is  the  variance  of  X(n)  which  is  assumed  to  be  the  same  for  all  n.  If, 
in  addition,  the  disorder  is  small,  we  obtain  from  Eq. (29) . 


2  2 
a  o 


Y  = 


8S2vsin20 


(45) 


For  other  cases  for  which  the  use  of  Eq . (27 ) ,  Eq.  (32)  or  Eq.(35)  is  more 
appropriate,  numerical  integrations  are  generally  required  using  an  assumed 
probability  density  for  X(n). 

The  computed  localization  factors  in  the  first  four  wave-passage  frequency 
bands  are  shown  in  Fig. 5(a)  and  (b)  for  6=1,  o=0.1  and  for  6=0.1,  o=0.1, 
respectively.  Except  for  frequencies  very  close  to  the  lower  boundary  of  a 
wave-passage  frequency  band  (not  shown  in  Fig. 5),  all  the  four  equations  (27), 
(29),  (32)  and  (35)  yield  essentially  the  same  results,  and  in  excellent 
agreement  with  those  obtained  from  Monte  Carlo  simulations,  also  shown  in  the 
figure.  That  the  values  of  the  localization  factors  are  below  0.02  is  an 
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Indication  of  weak  localization  in  these  cases.  At  the  lower  boundary  of  a 
wave-passage  frequency  band,  the  deformations  of  two  neighboring  spans  tend  to 
be  out  of  phase  (exactly  out  of  phase  for  an  ideal  periodic  structure)  and  the 
torsional  spring  at  the  supports  influence  the  beam  motion  to  a  great  extent. 
Therefore  even  a  small  disorder  of  the  torsional  springs  may  give  rise  to  strong 
localization.  On  the  other  hand,  the  deformations  of  two  neighboring  spans  are 
nearly  in  phase  at  the  upper  boundary  of  a  wave-passage  frequency  and  the 
torsional  springs  have  little  effect  on  the  beam  motion.  Hence,  their  disorder 
is  relatively  unimportant. 

The  Monte-Carlo  simulations  referred  to  above  were  carried  out  by 
multiplying  a  large  number  of  random  wave  transfer  matrices  at  a  given  frequency 
to  obtain  one  realization  of  the  wave  transfer  matrix  for  the  overall  N-span 
system  according  to  Eq . ( 17 )  and  thus  one  realization  of  tr(N,l).  The  values  of 
N_1£n  |  tr(N, 1)  J  for  different  realizations  were  averaged  over  a  sufficient 
number  of  realizations  to  obtain  an  estimate  for  the  localization  factor.  This 
procedure  is  different  from  the  one  adopted  by  Kissel16  in  which  an  average  is 
taken  of  a  large  number  of  £n|  t(n)  |  where  each  t(n)  is  the  transmission 
coefficient  of  one  realization  of  a  single  cell.  It  is  felt  that  our  simulation 
procedure  is  physically  more  meaningful.  In  our  simulations,  300  cells  were 
used  in  each  realization,  and  the  average  was  taken  over  2000  realizations. 

Fig. 6(a)  and  (b)  depict  the  computed  localization  factors  for  systems  with 
6=0.01,  o=0.25  and  6=0.01,  o=0.50,  respectively,  and  with  frequencies  falling 
within  the  lower  part  of  the  first  wave-passage  frequency  band.  It  is  seen  that 
the  closer  a  frequency  is  to  the  lower  boundary  of  the  wave-passage  frequency 
band,  the  more  Eqs.  (27)  and  (29)  overestimate  the  localization  factor,  because 
the  reflection  coefficients  become  larger  and  the  effect  of  multiple  reflection 
is  no  longer  negligible.  By  taking  the  reflections  in  only  one  cell  on  the 
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immediate  left  into  account,  Eq. (32)  yields  much  more  accurate  results  than 
those  of  Eqs.(27)  and  (29)  when  compared  with  simulation  results.  Fig.  6  shows 
that  for  a  beam  system  with  fixed  $  and  o,  localization  is  strong  at  a  frequency 
very  near  to  the  lower  boundary  of  a  wave-passage  frequency  band,  and  it  becomes 
weaker  as  the  frequency  increases.  This  implies  that  the  disorder  to  internal 
coupling  ratio  (o/S)  is  not  the  only  factor  in  determining  the  intensity  of 
localization,  as  once  suggested;  the  frequency  location  relative  to  a  wave- 
passage  band  is  at  least  another  factor.  Other  factors  might  be  also  important 
in  more  complicated  systems. 

Results  obtained  from  Eq.(27)  are  generally  more  accurate  than  those 
obtained  from  Eq.(29).  However,  exceptions  can  occur,  as  shown  in  the  higher 
frequency  region  in  Figs. 6(a)  and  6(b).  This  is  due  to  the  fact  that  both 
equations  are  approximate.  Although  Eq.(27)  essentially  retains  more  terms  in 
the  approximatation,  it  does  not  follow  that  those  additional  terms  will 
increase  the  accuracy  consistently. 

As  pointed  out  above,  localization  may  be  strong  at  frequencies  near  the 
lower  boundaries  of  wave-passage  frequency  bands,  even  when  the  disorder  is 
relatively  small.  Three  such  cases:  (a)  S=0.1,  w=232.5  rad/s,  (b)  5=0.01, 
w=318.7  rad/s  and  (c)  6=0.001,  w=339.72  rad/s  have  been  examined  in  more  detail. 
The  results  obtained  are  shown  in  Fig.  7 (a)-(c) ,  respectively,  with  varying 
standard  deviation  of  disorder  o.  In  all  these  cases  Eq . ( 32 )  is  clearly 
inadequate.  Eq. (35)  leads  to  satisfactory  results  for  cases  (b)  and  (c).  For 
case  (a)  in  which  the  disorder  is  also  greater,  Eq.(35)  is  still  inadequate 
although  it  is  much  superior  to  Eq. (32) .  The  use  of  even  higher  order 
approximations  is  suggested. 


Concluding  Remarks 

The  successive  improvement  in  accuracy  and  the  inclusion  of  structural 
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damping  are  two  important  features  in  the  proposed  perturbation  scheme  for  the 
computation  of  localization  factors  for  disordered  periodic  structures.  The 
first  feature  permits  a  unified  approach  to  any  disorder  to  internal  coupling 
ratio,  and  at  any  frequency  location  relative  to  the  lower  boundaries  of  the 
wave-passage  frequency  band.  They  are  shown  to  be  the  governing  factors  for  the 
degree  of  localization  in  the  illustrative  examples.  The  second  feature  permits 
the  consideration  of  practical  structures.  The  formulation  based  on  a  generic 
disordered  periodic  structure  is  also  an  advantage,  since  it  is  not  restricted 
to  a  specific  type  of  system,  and  it  facilitates  future  extension  to  the  case  of 
multiple  coupling  between  neighboring  cells. 
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Fig.  1  State  vector  representation;  (a)  Original  system,  (b)  Equivalent 
system. 
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Fig.  3  Schematic  representation  of  transmitted  waves  through  the  Interface 


Localization  factors  at  frequencies  In  wave-passage  frequency  bands 


Localization  factors  at  frequencies  In  the  lower  half  of  the  first 
wave-passage  frequency  band;  5  ■  0.01,  o  «  0.25 


Localization  factors  at  frequencies  In  the  lower  half  of  the  first 
wave-passage  frequency  band;  S  «  0.01.  o  -  0  5 


7(a)  Localization  factors  at  frequencies  near  the  lower  boundaries  of  the 
first  wave-passage  frequency  band;  5  -  0.1,  w  ■  232.5  rad/s 
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Fig.  7(c)  Localization  factors  at  frequencies  near  the  lower  boundaries  of  the 
first  wave-passage  frequency  band;  6  •  0.001,  w  -  339.72  rad/s 


EXACT  AND  APPROXIMATE  SOLUTION  TECHNIQUES  IN  NON-LINEAR  RANDOM  VIBRATION 

Consider  a  dynamical  system  governed  by 

Sft  Xj  =  Tj (X,t)  +  gjk(X,t)Fk(t); 

j  =  1,2, _ ,N;  k  =  1,2,... ,M  (1) 

in  which  t  is  the  time  variable,  Xj  are  components  of  the  system  response  vector 

X,  Fk(t)  are  random  processes  representing  the  random  excitations,  and  a 

repeated  subscript  in  a  product  indicates  a  summation.  Functions  fj  and  gjk 

are  generally  nonlinear;  however,  their  functional  forms  are  deterministic. 

Two  types  of  random  excitation  are  possible.  A  particular  Fk(t)  is  a 
multiplicative  (or  parametric)  excitation  if  its  accompanying  factor  9 j k ( X , t ) 
depends  explicitly  on  the  response  X;  it  is  an  additive  (or  external)  excitation 
if  its  accompanying  factor  does  not  depend  on  X.  In  other  words,  a 

multiplicative  excitation  appears  as  the  coefficient  of  an  unknown  in  the 
governing  equations  whereas  an  additive  excitation  appears  as  an  inhomogeneous 
term.  The  roles  played  by  the  two  types  of  random  excitations  are  fundamentally 
different. 

Exact  Solutions 

Exact  solutions  for  randomly  excited  nonlinear  systems  are  difficult  to 
obtain.  The  possibility  for  exact  solutions  exists  if  the  random  excitations 
can  be  idealized  as  Gaussian  white  noises.  In  such  a  case  the  response  of  a 
system  to  random  excitations,  represented  as  a  vector  in  a  state  space,  is  a 
Markov  vector.  The  probability  density  of  a  Markov  vector,  conditional  on  the 
knowledge  of  the  initial  state,  is  governed  by  a  parabolic  partial  differential 
equation,  called  the  Fokker-Planck  equation,  which  can  be  derived  from  the 
equations  of  motion.  Still,  exact  solutions  for  the  Fokker-Planck  equation  is 
known  only  for  very  special  first  order  nonlinear  systems,  which  are  of  little 
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interest  to  structural  engineers.  For  higher  order  systems,  solutions  have  been 
obtained  for  the  reduced  Fokker-Planck  equations;  namely,  Fokker-Planck 
equations  without  the  time  derivative  term.  The  unknown  in  a  reduced 
Fokker-Planck  equation  is  the  unconditional  probability  density  of  the  system 
response  when  it  reaches  the  state  of  statistical  stationarity.  This  can  occur 
if  the  fj  and  gj|<  functions  in  equation  (1)  do  not  depend  explicitly  on  t,  and 
if  energy  can  be  dissipated  from  the  system  to  balance  the  input  energy  from  the 
excitations. 

The  Fokker-Planck  equation  has  the  general  form  of 

ft  q  +  5xJ  GJ  =  0  (2) 

where  q  is  an  abbreviated  symbol  for  the  unknown  conditional  probability  density 
qx(x,t|  x0,t„),  and  Gj  is  the  probability  flow  in  the  jth  direction  given  by 

Gj  =  ~  \  (3) 

In  (2)  and  (3),  a  lower  case  symbol  represents  the  possible  value  of  a  random 
quantity.  For  example,  x  is  a  possible  value  of  X,  xj  a  possible  value  of  X j , 
etc.  This  convention  will  be  followed  throughout  the  paper.  In  equation  (3) 
Aj  and  Bj|<  are  called  the  first  and  second  derivate  moments  (also  known  as  the 
drift  and  diffusion  coefficients),  respectively,  which  can  be  obtained  from  the 
equations  of  motion  as  follows  using  the  well-known  rules  of  Wong  and  Zakai  [1]: 


Aj(x)  =  f j (x)  +  ttK£s  grs  ^  gj£(x) 

(4) 

Bjk(x)  =  2ttK£s  9j SL  9ks 

(5) 

In  (4)  and  (5)  each  <£S  is  a  constant,  representing  the  joint  spectral  density 
of  F$(t)  and  Fs(t),  which  are  now  restricted  to  be  Gaussian  white  noises  so  that 
an  exact  probabilistic  solution  can  be  obtained.  K^s  is  related  to  the  cross¬ 
correlation  function  of  F$(t)  and  -(t)  as 
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(6) 


E[FA(t1)Fs(tl)]  =  2ttK£s  S(t,  -  t2) 

An  exact  solution,  if  obtainable,  must  satisfy  the  initial  condition 

,  N 

Qx(x,t0|  x0,t)  =jniS(Xj  -  xj>0)  (7) 

and  appropriate  boundary  conditions.  If  the  total  probability  measure  must  be 
preserved,  then  the  boundary  conditions  must  be  natural,  namely, 


Gj  =  0  at  the  boundaries 


(8) 


We  shall  be  interested  in  the  stationary  state  solutions,  governed  by  the 
reduced  Fokker-Planck  equation 


Gj  -  Aj(x)p,(x)  -  jsfr  [BjkMM*)] 


(9) 


(10) 


where  ps(x)  is  the  stationary  probability  density  of  X  independent  of  t. 
Equation  (9)  is  solved  with  the  boundary  conditions  (8). 


1.  The  Case  of  Stationary  Potential 

The  special  case  of  stationary  potential  considered  by  Stratonovich  [2]  is 
one  in  which  the  probability  flow  vanishes  everywhere  not  only  at  the 
boundaries,  namely 

Aj(x)ps(x)  -  j  [8jk(x)Ps(x)  ]  =  0  (11) 
Without  loss  of  generality,  write  the  solution  for  (11)  as 

ps(x)  =  Cexp(-0)  (12) 
This  solution  form  guarantees  that  the  probability  density,  if  it  exists,  is 
non-negative  for  a  positive  normalization  constant  C.  Substituting  (12)  into 
(11),  one  obtains 

BJk  =  Bjk  ‘  2AJ  (13) 
These  are  N  equations  for  0.  If  a  consistent  0  function  can  be  obtained  from 
(13),  then  the  stochastic  system  belongs  to  the  class  of  stationary  potential. 
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For  a  first  order  system,  for  which  the  state  space  is  one-dimensional, 
namely  N  =  1,  the  probability  flow  must  vanish  everywhere  if  it  vanishes  at  the 
boundaries.  Then  (13)  reduces  to 


B  at  ’  ax  B  -  2A 


(14) 


Solving  for  $ 

♦  «  4nB  -  2  J  AB'1  dx  (15) 

Unfortunately,  as  to  be  explained  later,  oscillatory  dynamic  systems  do  not 
belong  to  the  class  of  stationary  potential.  Since  structural  systems  are 
usually  oscillatory  a  more  general  procedure  is  needed  and  will  be  discussed  in 
the  next  section. 


2  The  Case  of  Generalized  Stationary  Potential 

Circulatory  probability  flows  are  generally  present  in  a  multi-dimensional 
state  space.  To  provide  for  additional  flexibility  in  obtaining  a  solution  to 
equation  (9),  we  split  the  drift  and  diffusion  coefficients  as  follows  [3,4] 


Aj  .  +  Aj2) 

(16) 

BJk  •  »«>  +  B(j’ 

(17) 

Equation  (17)  retains  the  symmetry  condition  for  diffusion  coefficients;  i.e. 
Bjk  =  Bkj*  Substituting  (16)  and  (17)  into  (9),  one  obtains 


3 


Aj 


1) 


Ps 


*•»> 


Ps 


Let  Aj ( 1 )  be  associated  with  the  vanishing 
associated  with  the  circulatory  probability  flow, 
form  of  (12),  it  can  be  shown  that 


5JJ  [a)2)Ps]  *  0  (IS) 

probability  flow,  and  Aj(^) 
Then  for  ps  expressed  in  the 


jSj  *S2>  ■  *i2)  % 
Aj”  ■  A 


(19) 

(20) 
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3$ 


The  superscript  (j)  In  (20)  does  not  imply  a  summation.  There  are  a  total  of  n 
+  1  equations  fur  $  In  the  set  of  (19)  and  (20).  The  problem  is  solvable  if  a 
consistent  ♦  function  can  be  obtained  from  these  equations.  Such  a  solvable 
problem  is  said  to  belong  to  the  case  of  generalized  stationary  potential  [3,4]. 
It  is  of  interest  to  note  that  if  the  system  is  in  a  state  known  as  detailed 
balance  [e.g.  5,6],  then  the  two  parts  of  each  drift  coefficient  Aj(*)  and 
Aj (2)  correspond  to  the  so-called  Irreversible  and  reversible  parts, 
respecti vely.  However,  (19)  and  (20)  are  more  general,  and  are  applicable  to 
systems  not  in  detailed  balance. 

3  Applications 

The  above  solution  technique  will  now  be  applied  to  single-degree-of- 
freedom  second  order  systems  of  the  type 

2  +  iKZ.J)  «  0|t(2,i)Wk(t)  (21) 

where  ( t)  are  Gaussian  white  noises  with  delta  type  correlations;  i.e. 

E[Wk(t)W£(t+i)]  =  2nKla  6(1)  (22) 

Equation  (21)  is  equivalent  to  two  first  order  equations  in  the  form  of  (1) 

3fX,  =  X2  (23) 

^x2  =  -  *(Xi,X2)  +  0|c(X1,X2)  Wk(t) 

It  follows  from  (4)  and  (5)  that 

Ax  =  x2, 

A2  =  -  ^(Xx.Xj)  +  TTKk£  0£(X1,X2)g|-0k(X1,X2) 

011  =  012  =  011  =  02i  =  2lTKk£0k(X1  ,x2  )0£(xx  ,x2 ) 

Split  the  drift  and  diffusion  coefficients  into 

a!1)  =  0,  a!2)  =  x2 

^2  =  -  f(Xa,XS)  +  TlKk£  0£  8k  +  Xx/Xy 
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Ay/At 


B(2)  - 

°2  2 


b2 


2  » 


BU) 

D1  2 


-  - 


2  X 


1  * 

where  A  Is  an  arbitrary  function  of  x  =  xt  and  y  =  jx,,  and  Is  twice 
differentiable  with  respect  to  its  arguments.  The  physical  meaning  of  A  will  be 
discussed  later.  Then  application  of  (19)  and  (20)  results  in 


'  x*  Mr  +  <*x/Ay)  577  ‘  577  (Ax/Xy}  =  0 


3$ 


(24) 


£  •  £  B<‘> 

5x7  +  "KkA  0k  84  5x7  =  5x7  +  'l'(Xl,X2) 


(25) 


+  0k  ^  0a  -  (26) 

The  general  solution  of  (24)  is 

$(Xi  )X2)  =  -  An  Ay  +  4>0(A)  (27) 

where  <^0 (A)  is  an  arbitrary  function  of  A.  In  order  that  (25)  and  (26)  can  also 
be  satisfied,  the  and  8k  functions  must  be  related  as 


'HXi  ,X2)  =  1TX2K|<£  0k  0£  [Ay  <Do  1  (A)  -  Ayy/Ay] 

-  I*kt  Oa  g|7  ak  +  Xx/Ay  .  D(x,)  e*t(X)/Ay  (28) 

where  D(xJ  is  an  arbitrary  function  of  x2  alone.  Equation  (28)  provides  a 
restriction  between  the  constitutive  law  of  the  system,  as  represented  by  \Ji  and 
0k  on  one  hand,  and  the  spectral  densities  Kk£  on  the  other,  so  that  the  exact 
solution  for  the  stationary  probability  density  is  obtainable  using  the  present 
procedure. 

The  well-known  case  [7,8]  of  linear  damping,  nonlinear  spring,  and  an 
additive  Gaussian  white  noise  excitation 

Z  +  82  +  g(Z)  =  W(t),  8  >  0  (29) 

corresponds  to  a  single  0  =  1  and 

’KXi.Xj)  =  0x2  +  g(x2)  (30) 
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(31) 


Then  letting  D(xJ  =  0  and 

x  1  i  *i 

A  =  y  +  /  g(u)du  =  7x,  +  J  g(u)du 
0  ^  0 

equation  (28)  reduces  to 

MX)  =  |u  (32) 

where  K  is  the  spectrum  of  the  white  noise.  Furthermore,  $(A)  =  $„(A)  in  this 
case.  Thus 

ps(x1,x2)  =  C  exp(-  ^  A) 

=  C  exp{-  ^  [|  x|  +  f1  g(u)du ]}  (33) 

which  is  a  valid  probability  density  if  A  Is  positive  for  large  |  xx  |  . 

Without  specifying  the  form  for  the  A  function,  but  still  letting  D(xx)  =  0 
and  restricting  to  a  single  additive  excitation,  one  obtains  from  (28) 

<j,0'(X)  s  Mxi>xi.?  .t,  .ftx^  (34) 

Note  that  the  right  hand  side  of  (34)  must  be  a  function  of  A  alone.  If  this  is 
satisfied,  then  (34)  can  be  used  to  obtain  the  required  <MA). 

The  solution  obtained  by  Caughey  and  Ma  [9]  corresponds  to 

Y(Xi,x2)  a  Ax/Ay  +  ttKXj  [Ayh(A)  -  Ayy/Ay]  (35) 

where  h(  )  is  an  arbitrary  function  subject  to  a  rather  general  restriction  to 
be  specified  later.  Substitution  of  (35)  into  (34)  results  in 

$o'(A)=h(A)  (36) 

and  a  stationary  probability  density 

A 

PstXi.Xi)  =  CA  exp  [-  /  h(v)dv]  (37) 

y  o 

The  A  and  h  functions  must  be  such  that  the  right  hand  side  of  (37)  is 
normalizable.  This  extremely  simple  result  is  obtained  under  a  rather 
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artificial  assumption  that  the  constitutive  law  of  the  system  itself,  namely 
'l'(Xi.Xj)  is  dependent  on  the  spectrum  K  of  the  excitation. 

The  physical  meaning  of  the  A  function  can  be  revealed  by  considering  the 
deterministic  case  of  undamped  free  motion  governed  by 

x  +  Ax/Ay  =  0  (38) 

Since  x  =  (1/2) (dx*/dx)  =  dy/dx,  the  above  equation  can  be  replaced  by 

dA  =  Axdx  +  Aydy  =  0  (39) 

Therefore,  A(x,y)  is  not  truly  arbitrary;  it  represents  the  total  energy  in  the 
conservative  system  (38).  It  should  be  noted,  however,  that  this  physical 
interpretation  must  be  modified  if  multiplicative  random  excitation  appears  in 
the  coefficient  of  the  first  derivative  term  in  a  second  order  differential 
equation. 

Equations  (19)  and  (20)  have  been  applied  to  obtain  new  solutions  for 
systems  under  both  additive  and  multlpl icative  excitations  [4,10].  The 
following  two  examples  will  illustrate  some  of  the  possibilities. 

The  first  example  is  described  by 

Z  +  [a  +  j  (2  +  w0Z)  +  W2 (t) ]Z 

+  Wo  [1  +  Wx(t) JZ  =  W,(t)  (40) 

where  Wx ( t) ,  W2(t)  and  W3(t)  are  independent  Gaussian  white  noises  with  spectra 
<ii,  <22  and  K3j,  respectively.  Under  the  restriction  w0‘>K11=K22  the 
stationary  probability  density  for  Z=X2  and  Z=X2  is  given  by  [14] 

Ps(*i,x2)  s  C(A  +  ^  exp(-  A)  (41) 

where  y  -  ^  (ttK2 2  ^  ~  i  <jj)  and  A  =  ^  x2  +  ^  Wo  x2 

TT<2  2 

Equation  (41)  is  a  generalization  of  a  result  obtained  previously  by  Dlmentberg 
[10]  for  the  special  case  K33  =  0,  namely,  when  the  additive  random  excitation 
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Is  absent,  In  which  case,  the  value  for  y  is  restricted  to  y  <  1;  otherwise,  the 
right  hand  side  of  (41)  is  not  normalizable.  This  can  occur,  of  course,  only  if 
a  <  0  and  k22  >  |  a  |  /it.  It  means  that  the  trivial  solution  Z=2=0  is  unstable  in 
probability  when  k22  >  |  a  |  it. 

The  second  example  is  an  Interesting  nonlinear  oscillator  capable  of 
exhibiting  limit  cycle  and  bifurcation  behaviors,  governed  by  [11]. 

2  +  z2[0  +  -T^TlJ  z  +  [1  +  W(t)]  Z  =  0  (42) 

z  +z 

where  0  >  0,  and  W(t)  is  a  Gaussian  white  noise  multiplicative  excitation  with  a 
spectral  density  K.  Without  the  random  excitation,  the  equation  has  a  non¬ 
trivial  solution,  corresponding  to  a  stable  limit-cycle  if  a  <  0.  If  a  >  0, 
then  the  only  stable  configuration  is  the  trivial  Z  =  0. 

The  stationary  probability  density  for  Z=X2  and  ?=X2  is  found  to  be 

Ps(Xi,x2)  *  C(x*  +  x!)(2a/TTK)  exp[-  ^(xj  +  x!)]  (43) 

Equation  (43)  implies  that  if  a  >  0  then  a  nontrivial  stationary  probability 
density  exists  only  if  K  >  2a/iT.  Thus,  2a/ir  represents  a  critical  spectral 
level  for  the  multiplicative  excitation  at  which  bifurcation  occurs,  and  below 
which  the  trivial  solution  is  stable  in  probability.  When  a  stationary 

probability  density  of  the  system  response  exists  (K  >  2a/ir) ,  its  shape  depends 
on  whether  a  2  0  or  a  <  0.  It  can  be  shown  that  if  a  i  0,  then  the  probability 

density  ps  has  a  peak  at  x2  =  x2  =  0.  On  the  other  hand,  if  a  <  0,  then  ps  has 

2  2  2  2 
a  maximum  at  x,  +  x2  =  4  |  a  |  /0  and  it  tends  to  zero  as  xl  +  x2  approaches  to 

either  zero  or  infinity.  This  Implies  that  the  probability  density  has  a  shape 

shown  In  Fig.  1,  and  the  most  probable  values  for  the  response  lie  on  the  circle 

Xi  +  x2  *  4  |  a  |  /0,  which  is  the  stochastic  analogue  of  the  limit  cycle. 

To  the  knowledge  of  the  writer,  the  class  of  generalized  stationary 
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potential,  satisfying  (19)  and  (20),  includes  all  the  stochastic  systems  whose 
exact  probability  solution  in  the  stationary  state  are  known  to  date.  When  both 
multiplicative  and  additive  random  excitations  are  present,  restriction  on  the 
system  parameters  and  excitation  spectral  levels  are  generally  required  to 
obtain  such  solutions. 

Approximate  Solutions 

When  both  multiplicative  and  additive  random  excitations  are  present, 
restrictions  usually  must  be  imposed  on  the  relationship  between  the  noise 
spectra  and  some  system  parameters  in  order  to  obtain  an  exact  stationary-state 
solution  for  the  response  of  a  nonlinear  system  at  the  present  time.  In 
practical  engineering  applications,  the  excitation  spectral  levels  cannot  be  so 
restricted;  therefore,  approximate  procedures  are  often  needed  to  obtain  the 
solutions. 

1  Cumulant-Neglect  Closure 

One  versatile  scheme  to  obtain  approximate  solutions  is  the  cumulant- 
neglect  closure  [14,15]  which  can  be  applied  easily  if  random  excitations  are 
Gaussian  white  noises  or  filtered  Gaussian  white  noises.  Return  to  (1)  and 
convert  these  equations  to  the  I  to  type  stochastic  differential  equations  as 
follows  [e.g.  16] 


dXj(t)  =  mj(X,t)dt  +  oj£(X,t)dB£(t)  (44) 

where  B£( t )  are  independent  unit  Wiener  processes,  and 

mj(X,t)  =  Aj (X , t )  (45) 

Oj*(X,t)  okt(X,t)  =  Bjk(X,t)  (46) 

The  right  hand  sides  of  (45)  and  (46)  are  the  same  drift  and  diffusion 
coefficients  in  (4)  and  (5),  but  with  x  replaced  by  X.  However,  these 
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coefficients  can  now  depend  on  t,  and  the  results  in  this  section  are  not 
restricted  to  the  stationary  response.  One  can  then  obtain  equations  for  the 
statistical  moments  for  X(t)  by  using  Ito's  differential  rule  [17] 


—  E[F(X,t)]  =  E[—  +  (mj  5E_  +  I  o j £ok£  -^— )]  (47) 

dt  at  8Xj  2  8Xj9Xk 

and  by  letting  F  be  X^ ,  X^Xr,  XiXrXs,  etc.  For  nonlinear  systems,  the  equations 
for  statistical  moments  form  an  infinite  hierarchy,  and  they  can  only  be  solved 
approximately  using  a  suitable  closure  scheme. 

The  statistical  moments  are  related  to  cumulants  as  follows  [2,18]: 

E[Xj]  =  Kt[Xj] 

E[XjXk]  =  <2[Xj,Xk]  +  ^[Xj^pCk] 

E[XjXkX£]  .  K,[Xj,Xk,Xt]  +  3  {  K  2 [X j ]k2 [Xk,X£]}s 
+  x1[Xj]x1[Xk]  k2[X£] 

ECXjX^Xfj,]  =  K^Xj.Xu.Xa.Xh,]  +  3{KI[Xj,Xk]K2[X£,Xm]}s 
+  4{K1[Xj]K2[Xk,X£,Xfn]}s  +  6{K1[Xj]K1[Xk]x2[X£,Xm]}s 
+  X  i  [X  j  ]k  2  [Xk]x  t  [X£]k  2  [X,n] 

.  (48) 

where  Xj[  ]  denotes  a  jth  cumulant,  and  {  }s  indicates  a  symmetrizing  operation 
with  respect  to  its  arguments;  namely,  taking  the  arithmetic  mean  of  different 
permuted  terms  similar  to  the  one  within  the  braces.  In  a  cumulant-neglect 
closure,  all  the  cumulants  higher  than  a  given  order  N  are  set  to  zero.  The 
statistical  moments  of  an  order  higher  than  N  can  then  be  expressed  in  terms  of 
lower  order  moments  using  (48).  Very  accurate  results  have  been  obtained  using 
this  scheme  [14]  as  long  as  the  spectral  levels  of  multiplicative  random 
excitations  are  not  too  high  so  that  the  system  is  far  from  stochastic 
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bifurcation  or  instability  [19,20]. 

It  Is  of  Interest  to  note  that  for  a  random  process  defined  on  (--,  -) , 
neglecting  xn  for  n  >  2  Is  equivalent  to  the  so-called  Gaussian  closure;  namely, 
approximating  this  random  process  by  a  Gaussian  process.  Gaussian  closure, 
however,  is  the  simplest  and  of  the  lowest  order  among  non-trlvial  cumulant- 
neglect  closures.  Furthermore,  if  the  excitations  are  white  noises  or  filtered 
white  noises,  then  Gaussian  closure  leads  to  the  same  results  as  statistical 
equivalent  linearization.  A  review  of  recent  contributions  to  the  equivalent 
linearization  schemes  has  been  given  by  Spanos  [21].  Another  closure  scheme, 
using  a  truncated  Gram-Charl ier  or  Edgeworth  series  [22,23]  has  a  similar 
accuracy  as  that  of  cumulant  closure;  however,  it  is  more  difficult  to  apply  to 
multi -degree-of-freedom  systems  when  a  multidimensional  Gram-Charl ier  series  is 
required. 

2  Generalized  Equivalent  Linearization 

The  development  of  generalized  equivalent  linearization  scheme  by  Bruckner 
and  Lin  [24]  was  motivated  by  two  reasons:  (1)  a  cumulant-neglect  closure  of  an 
order  N>2  does  not  guarantee  a  non-negative  probability  density  everywhere,  (2) 
the  traditional  equivalent  linearization,  which  leads  to  the  same  result  as  that 
of  Gaussian  closure,  replaces  the  original  nonlinear  system  with  a  linear  system 
devoid  of  parametric  excitation;  therefore,  it  is  unsuitable  to  treat  the  case 
when  random  parametric  excitations  are  also  present. 

In  a  generalized  equivalent  linearization  scheme,  the  original  Ito  system, 
(44)  is  replaced  by  a  linear  Ito  system  as  follows: 

dXj(t)  =  (ajkXk  +  bjjdt  +  (cjrsXr  +  djs)dBs(t)  (49) 
The  equations  for  the  first  and  second  order  statistical  moments  are  closed; 
they  may  be  determined  from 
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5t  ECxjJ  -  ajk  E[Xk]  +  bj 


(50) 


ECxjxr]  s  ars  5[XjXs]  +  ajk  E[XrXk] 

+  cjk4  crp£  E[XkXp]  +  br  E[Xj]  +  bj  E[Xr] 

+  Cjk£  df£  E[Xk]  +  crp£  dj£  E[Xp]  +  dj£  dr£  (51) 

The  constants  ajk,  bj,  Cjrs  and  djS  are  determined  by  minimizing  the  mean-square 

errors 

Cj  =  E{[mj(X)  -  ajkXk  -  bj]*}  (52) 

£jr  =  5{[Xjfr>r(X)  +  Xpmj(X)  +  Oj£(X)  Op£(X)  -  a|~sXjX5 
-  ajkXrXk  -  Cjk£  c r p £  XkXp  -  brXj  -  bjXr 

■  cjk4  dr£  xk  "  crp i  dj SL  xp  ~  dj£  dr£J  }  (53) 

To  carry  out  the  minimization  process,  the  evaluation  of  certain  moments  higher 
than  the  second  order  Is  still  required.  However,  when  the  optimal  constants  so 
obtained  are  substituted  into  Eq.(49),  the  linearized  system  does  possess  a  non¬ 
negative  probability  density. 

3  Stochastic  Averaging 

Method  of  stochastic  averaging  was  developed  by  Stratonovich  [2]  based 
primarily  on  physical  arguments.  It  was  later  justified  and  interpreted 
rigorously  by  Khasminskii  [25],  Papanicolaous  and  Kohler  [26].  The  first 
application  to  structural  systems  was  due  to  Arlaratnam  in  the  anlaysis  of 
column  stability  under  random  axial  load  [27].  Lin  and  coworkers  have  used  the 
procedure  to  Investigate  the  effects  of  atmospheric  turbulence  on  the  motions  of 
helicopter  blades  [e.g.  28,29,30,31],  and  of  long-span  bridges  [32,33]. 

Namachichivaya  and  Arlaratnam  [34]  and  Namachchivaya  [35,36]  applied  the  method 
to  Hopf  bifurcation  in  the  presence  of  stochastic  excitations.  A  review  of  the 
general  subject  and  many  other  contributions  was  given  by  Roberts  and  Spanos 
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[37]. 


According  to  Stratonovich  [2]  the  response  of  a  dynamic  system  to  random 
excitation  may  be  approximated  by  a  Markov  vector,  if  the  correlation  times  of 
all  the  excitations  are  short  compared  with  the  relaxation  time  of  the  system. 
Then  the  coefficients  in  the  corresponding  Ito  type  differential  equations  (44) 
can  be  obtained  from  the  original  physical  equations  (1)  as  follows: 

mj  =  9jk(X.t)  g£s(X»t  +t) 

E[Fk(t)  Fs(t  +  i)]di  (54) 

OD 

°j£  °ki  =  f  «jr(X,t)gjcs(X,t+T)  E[Fr(t)  Fs(t+x)]dt  (55) 

It  is  of  interest  to  note  that  the  right  hand  sides  of  (54)  and  (55)  reduce  to 
the  same  forms  as  those  of  (4)  and  (5),  respectively,  if  the  random  excitations 
are  delta-correlated. 

If  the  right  hand  sides  of  the  physical  equations  (1)  are  small  in  some 
sense,  implying  that  the  state  variables  are  slowly  varying,  then  the  results 
obtained  in  (54)  and  (55)  can  be  further  simplified  by  time-averaging  over  t. 
This  additional  time  averaging  not  only  makes  the  drift  and  diffusion 
coefficients  independent  of  time,  but  very  frequently  it  also  renders  a  group  of 
state  variables  to  become  uncoupled  from  the  others,  thus  permitting  a  reduction 
of  the  dimensions  of  the  state  space  being  considered.  The  latter  advantage  is 
very  significant  from  an  analytical  point  of  view  and  is  the  main  reason  for  the 
popularity  of  the  stochastic  averaging  procedure.  If  the  system  is  linear  or 
slightly  nonlinear,  the  unknowns  In  the  original  physical  equations  can  be 
transformed  to  polar-coordinate  variables  corresponding  to  slowly  varying 
amplitudes  and  phases,  making  the  time-averaging  applicable.  Alternatively,  the 
unknowns  can  be  transformed  to  complex  random  processes  in  conjugate  pairs.  The 
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absolute  magnitude  of  a  complex  random  process  is  equivalent  to  an  amplitude 
process,  and  the  ratio  between  the  imaginary  and  real  parts  is  the  arctangent  of 
a  phase  process.  This  alternative  scheme  was  used  by  Ariaratnam  and  Tam  [38]  to 
linear  systems  with  many  degrees  of  freedom.  Bruckner  and  Lin  [39]  have  shown 
that  it  can  also  be  applied  advantageously  to  nonlinear  cases,  and  when  used  in 
conjunction  with  the  cumulant-neglect  procedure,  the  number  of  equations  to  be 
solved  simultaneously  is  reduced  greatly. 

Since  time-averaging  is  often  the  key  to  make  a  problem  solvable,  it  is  of 
interest  to  comment  on  two  cases  for  which  time-averaging  are  still  applicable 
without  requiring  that  the  right  hand  sides  of  all  the  first  order  physical 
equations  be  small.  The  first  case  is  a  weakly  damped  single-degree-of-freedom 
oscillator  with  a  strongly  nonlinear  stiffness,  governed  by 

X  =  -  g (X)  +  ef (X , X)  +  e^  h1 (X,X)Fi (t)  (56) 

where  e  is  a  small  positive  number  and  for  simplicity  the  F ^ (t)  are  assumed  to 
be  white  noises.  Introduce  an  energy  process  [2] 


V  =  jX2  +  G(X)  =  ^X*  +  (  9(u)du  (57) 

Equation  (56)  can  be  replaced  by  the  following  two  first  order  equations 

X  =  ±V  2(V  -  G)  (58) 

V  =  ±  V  2(V-G)  {ef (X,±  /  2(V-G)  ) 

+  e^hi(X,±/TrV^)  )Fi(t)}  (59) 

These  two  equations  can  be  converted  to  the  following  I  to  type  equations 

dX  =  ±  V  2 ( V-G )  dt  (60) 

dV  =  emdt  +  o  dB(t)  (61) 

where  according  to  (4),  (5),  (45)  and  (46) 

m  =  ±  ✓  2(V-G)  {f  +  (ttK1  j)  hj  ^j[±  /  2(V-G)hi  ]]  (62) 
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o  =  V  4uKij(V-G)hihj  (63) 

The  unknown  V  in  (65)  is  slowly  varying,  but  the  unknown  X  in  (60)  is  not. 
Invoking  a  theorem  due  to  Khasminskii  [40]  Zhu  [41]  performed  a  time-averaging 


on  the  right-hand-side  of  (65)  as  follows: 

m  -  }  /{f  +  irKij  hj  [±  rm1 G)  h-j  ]}dX  (64) 

R 

o2=  ±  J  /  2tiK1i/2(V-G)hihj  dX  (65) 

1  R 

where  integration  on  t  has  been  replaced  by  integration  on  X  using  (60),  R  is 
the  integration  domain  defined  as  V  £  G,  and 

T  «  ±  /  [ 2 ( V-G ) ] "^  dX  (66) 

R 

Since  the  integrands  are  expressed  in  terms  of  ±/2(V-G) ,  the  value  of  each 
integral  is  obtained  by  adding  the  results  from  separate  integrations  for 

the  plus  and  the  minus  signs.  Substituting  m  and  o  for  m  and  o,  respectively, 

in  (61),  we  obtain  an  Ito  equation  for  V  which  is  devoid  of  X.  Similar  analyses 

have  also  been  carried  out  by  Spanos  [42]  and  Roberts  [43]  for  systems  under 

only  additive  random  excitations. 

Another  exceptional  case,  considered  by  Namachchi vaya  and  Lin  [44],  is 
described  by  the  following  equations  of  motion 

Zi  +  2a j  Zx  +  nUi  =  f i (Z,Z,y)  +  A, j (p)ZjW(t) 

Z2  +  2a2Z2  +  Q22z2  =  f2(Z,Z,p)  +  A2  j  (p)ZjW(t)  (67) 

where  p  is  a  parameter,  nj2=wj*(p)+  a^2(p),  2nd  are  damping  parameters  and 
natural  frequencies,  respectively,  and  fj  are  only  weakly  nonlinear.  As  p 
increases,  the  first  mode  becomes  close  to  instability  while  the  second  mode 
remains  highly  damped.  We  then  transform  Z2  and  to  the  amplitude  and  phase 
counterparts,  and  at  the  same  time  make  the  substitution 
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(68) 


Z2  =  [cosw,t  +  wj1a2sinw2t)Y1  +  (w2 1sinw2t)Y2]e  ajt 

Z2  =  [-(wlViIsliWttJYj  +  (cosw2t-w21a2sinw2t)Y2]e'rilt 
The  new  equations  satisfy  the  condition  for  time  averaging,  and  the  resulting 
Ito  equation  for  the  amplitude  of  the  first  mode  is  devoid  of  other  variables. 
This  method  was  developed  on  the  basis  of  a  theorem  due  to  Papanicolaou  and 
Kohler  [45].  It  can  be  applied  to  cases  Involving  a  very  large  number  of 
degrees  of  freedom.  However,  if  more  than  one  mode  can  become  critical,  then 
transformation  to  conjugate  pairs  of  complex  processes  is  simpler  [46]. 

4  Energy  Dissipation  Balancing  [47] 

Consider  a  single-degree-of-freedom  system,  governed  by  equations  of  the 
type  of  (23).  An  exact  stationary  probability  solution  can  be  found  if  the  i|i 
and  0  functions  are  related  according  to  (28).  Restricting  to  D(x2)  =  0  and  to 
A  of  the  form  of  (31),  relation  (28)  is  simplified  to 

<KXi,x2)  =  ™2Kk£0k0£$' (A)  -  +  9(xi)  (69) 

This  equation  specifies  a  class  of  solvable  problems. 

Given  a  problem  with  a  vji  function  not  in  the  form  of  (69),  it  is  possible 
to  write 

Y(Xi,X2)  =  U(X2,X2)  +  g(Xx)  +  TTKk£0£  gj-  0k  (70) 

Then,  using  (12)  as  an  approximate  solution,  the  error  involved  is 

e  =  u(xl(xs)  -  nx2  Kk£0k0£<t>  +  2uKk£  9 £t~  8k  (71) 

OX  j 

A  reasonable  criterion  for  the  choice  of  $'(A)  is  to  equate  the  average  energy 
dissipations  between  the  given  system  and  the  substituting  system  [47]  for  which 
the  exact  solution  is  obtainable;  namely  letting 

/  x2 €  PsU^xJdXidx,  =  0  (72) 

where  ps  Is  chosen  to  be  the  stationary  probability  density  of  the  substituting 
system. 
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As  an  example,  consider  a  van  der  Pol  type  oscillator 


Z  +  ai  +  0ZJ2  +  yV  +  Z  +  ZWj(t)  +  2w,(t)  +  W,(t)  (73) 

By  equating  Z  =  Xt  and  1  =  Z2  and  comparing  with  (23),  we  can  identify 

2  3 

that  0X  =  x»,  0,  =  x2,  0j  =  1  and  ^x^x,)  =  ax2  +  Bx2x2  +  yx2  +  x2 


9 

TTKkJ.9)i  9k  =  ^(^21*1  +  K22x2  +  K23)  (74) 

3x2 


From  (70), 


vKxltx2)  =  u(x2,x2)  +  g(x2 )  +  ir^iXi+K.jXj+Kja) 

Therefore, 

g(Xi)  =  x2(l  -  itK21)  -  ttK2 j 

2  3 

u(x2,x2)  =  (a  -  ttK22)x2  +  Bx2x2  +  yx2 
It  follows  from  (31)  that 


(75) 

(76) 


A  *  \*\  +  j(l  -  1K»i)(Xx  -  x1>0)2  (77) 

where  xlj0  =  irK2J/(l  -  nK21).  In  order  for  A  to  be  always  positive,  it  requires 
that  1  -  tiK12  >  0  which  is  a  necessary  condition  in  probability. 

The  integration  in  (72)  can  be  simplified  by  a  transformation  x2  = 
[2A/(1-ttK21)]^cos0  +  xl>0f  and  x2  =  (2A)^  sin0,  resulting  in 


0'(A)  =  Ha+vKuL+  (g  +  3y)A _ (78) 

2tt(  1  — TtK 2 2)  (K21x2  jO^KjiXj  >0+Kj2)  +h[Kj 2  +3K2  2 (  1-ttK2 2  ]A 

The  approximate  stationary  probability  denisty  is 

Ps(*i»0  =  c  exP  C-  /  0'  (A)dAJ  (79) 

The  reader  is  referred  to  [47]  for  details. 

It  is  of  interest  to  note  that  If  only  the  additive  random  excitations  are 

present  then  the  approximate  probability  density  obtained  on  the  basis  of  energy 

dissipation  balancing  coincides  with  that  obtained  from  a  least  mean-square 
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error  criterion  proposed  by  Caughey  [48];  however,  the  two  results  are  different 
if  multiplicative  random  excitations  are  present.  In  the  latter  case,  the 
energy  dissipation  criterion  leads  to  correct  bifurcation  and/or  instability 
conditions  in  some  known  examples  while  the  least  mean-square  error  criterion  or 
stochastic  averaging  [47,49]  do  not. 


Concluding  Remarks 

The  length  limitation  of  this  report  does  not  permit  a  thorough  review  of 
all  important  advances  which  took  place  in  recent  years  toward  solving  nonlinear 
random  vibration  problems.  It  is  therefore  appropriate  to  suggest  additional 
materials  for  interested  readers.  The  earlier  review  by  Caughey  [50]  remains 
most  authoriative  for  works  prior  to  1971.  A  large  portion  of  a  survey, 
entitled,  "Methods  of  Stochastic  Dynamics"  is  devoted  to  nonlinear  stochastic 
systems.  For  the  very  useful  analytical  tools  of  equivalent  linearization  and 
stochastic  averaging,  [21]  and  [37]  are  again  suggested. 
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density  (with  the  right- front  quarter  removed) 
analogue  of  limiting-cycle  behavior  of  a  non- 
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PROPOSED  FUTURE  RESEARCH 


Continuation  of  the  present  research  may  proceed  in  several  directions.  We 
have  investigated  the  localization  effects  due  to  random  variabilities  of 
material  properties  and  manufacturing  processes  for  disordered  structural  units 
with  mono-coupling,  which  permits  only  one  type  of  wave  to  propagate  from  cell 
to  cell  in  the  structure.  It  must  be  generalized  to  the  case  of  multi-coupling, 
because  the  truss-type  structural  networks  are  capable  of  transmitting  as  many 
as  24  types  of  structural  waves.  The  generalization  is  a  major  challenge,  since 
multi-wave  transmission  permits  re-composition  or  dispersion  among  different 
wave  modes.  The  computation  of  the  localization  factor  in  the  multi-wave  case 
is  analogous  to  the  computation  of  the  so-called  largest  Lyapunov  exponent  for  a 
dynamical  system  defined  in  a  large  dimensional  space,  which  is  still  an 
unsolved  problem.  In  order  to  attempt  a  break-through,  both  analytical  and 
simulation  approaches  may  have  to  be  employed. 

Although  we  have  developed  several  rather  general  schemes  for  obtaining 
exact  and  approximate  solutions  for  randomly  excited  nonlinear  structures.  The 
application  to  aerospace  structures  have,  by  and  large,  been  limited  to  single- 
degree-of-freedom  systems.  Recognizing  the  fact  that  idealized  conditions  for 
exactly  solvable  systems  can  seldom  be  met  for  multi-degree-of-freedom  aerospace 
structures,  approximate  solutions  are  the  norm  rather  than  exceptions.  However, 
the  class  of  generalized  stationary  potential  is  too  wide  when  applied  to  the 
case  of  multi-degree-of-freedom  systems;  and  the  accuracy  of  an  approximate 
solution  for  this  case  on  the  basis  of  dissipation  energy  balancing  may  be 
questionable.  It  is,  therefore,  necessary  to  reduce  the  size  of  the  class  in 
order  to  obtain  meaningful  approximations. 

Finally,  a  new  stochastic  structural  control  algorithm  should  be  developed 
on  the  basis  of  controlling  wave  motions  in  disordered  structures.  This 
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approach  is  different  from  what  has  usually  been  used  in  the  past  in  which 
emphasis  is  placed  on  controlling  the  response  in  various  normal  modes.  We 
believe  that  the  modal  approach  is  inappropriate  for  large  space  structures 
exposed  to  essentially  transient  excitations,  for  which  there  is  not  sufficient 
time  for  the  modal  response  to  be  generated  fully.  A  wave  propagation  approach 
is  more  direct  and  more  effective. 
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